
==== Front
Genome Biol Evol
Genome Biol Evol
gbe
Genome Biology and Evolution
1759-6653
Oxford University Press UK

39235041
10.1093/gbe/evae194
evae194
Article
AcademicSubjects/SCI01130
AcademicSubjects/SCI01140
Assessing Mechanisms of Potential Local Adaptation Through a Seascape Genomic Approach in a Marine Gastropod, Littoraria flava
https://orcid.org/0000-0003-2312-5087
Cortez Thainá Departamento de Genética e Biologia Evolutiva, Instituto de Biociências, Universidade de São Paulo, São Paulo, Brazil

https://orcid.org/0000-0002-1854-0692
Sonoda Gabriel G Departamento de Genética e Biologia Evolutiva, Instituto de Biociências, Universidade de São Paulo, São Paulo, Brazil
Laboratório de Toxinologia Aplicada, Instituto Butantan, São Paulo, Brazil

https://orcid.org/0000-0003-4885-901X
Santos Camilla A Departamento de Genética e Biologia Evolutiva, Instituto de Biociências, Universidade de São Paulo, São Paulo, Brazil
Tree of Life, Wellcome Sanger Institute, Cambridge, CB10 1SA, UK

https://orcid.org/0000-0002-1302-5261
Andrade Sónia Cristina da Silva Departamento de Genética e Biologia Evolutiva, Instituto de Biociências, Universidade de São Paulo, São Paulo, Brazil

Yeh Shu-Dan Associate Editor
Thainá Cortez and Gabriel G Sonoda considered joint first authors.

Corresponding authors: E-mails: thainacortez@usp.br; soniacsandrade@ib.usp.br.
9 2024
05 9 2024
05 9 2024
16 9 evae19413 7 2024
20 9 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of Society for Molecular Biology and Evolution.
2024
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Understanding the combined effects of environmental heterogeneity and evolutionary processes on marine populations is a primary goal of seascape genomic approaches. Here, we utilized genomic approaches to identify local adaptation signatures in Littoraria flava, a widely distributed marine gastropod in the tropical West Atlantic population. We also performed molecular evolution analyses to investigate potential selective signals across the genome. After obtaining 6,298 and 16,137 single nucleotide polymorphisms derived from genotyping-by-sequencing and RNA sequencing, respectively, 69 from genotyping-by-sequencing (85 specimens) and four from RNA sequencing (40 specimens) candidate single nucleotide polymorphisms were selected and further evaluated. The correlation analyses support different evolutionary pressures over transcribed and non-transcribed regions. Thus, single nucleotide polymorphisms within transcribed regions could account for the genotypic and possibly phenotypic divergences in periwinkles. Our molecular evolution tests based on synonymous and non-synonymous ratio (kN/kS) showed that genotype divergences containing putative adaptive single nucleotide polymorphisms arose mainly from synonymous and/or UTR substitutions rather than polymorphic proteins. The distribution of genotypes across different localities seems to be influenced by marine currents, pH, and temperature variations, suggesting that these factors may impact the species dispersion. The combination of RNA sequencing and genotyping-by-sequencing derived datasets provides a deeper understanding of the molecular mechanisms underlying selective forces responses on distinct genomic regions and could guide further investigations on seascape genomics for non-model species.

natural selection
marine invertebrate
transcriptome
non-synonymous and synonymous polymorphisms
==== Body
pmcSignificance

The evolutionary mechanisms underlying the distribution of genetic variation and responses to changing environments in marine organisms remain incompletely understood, despite concerning indicators of climate changes. In this study, we employed a widely distributed gastropod species commonly found on tropical Western Atlantic rocky shores and mangroves to investigate this issue. By combining population genomics and transcriptomics data, we discovered distinct selective pressures acting on coding and noncoding regions. Moreover, genetic structuring between localities was observed, driven by significant contributions from both synonymous and noncoding sequence variants. Our research provides valuable insights into the intricate mechanisms governing the response of a non-model species to a continually changing environment.

Introduction

A more expansive understanding of the connection between genetic structure and local adaptation has become a primary goal in seascape genomic studies (Liggins et al. 2019). This objective can be achieved by combining genomic approaches with techniques capable of detecting associations between environmental factors and the spatial distribution of genomic diversity (Bay and Palumbi 2014; Segovia et al. 2020). By analyzing more loci across the genome, a higher number of regions potentially affected by selection can be accurately identified (Liggins et al. 2019).

Genomic approaches such as RNA sequencing (RNA-Seq) (Gleason and Burton 2015; Huang et al. 2016) and genotyping-by-sequencing (GBS) (Segovia et al. 2020) allow for the identification of thousands of single nucleotide polymorphisms (SNPs) without prior genetic information from the species (Gagnaire et al. 2012; Tepolt and Palumbi 2015). RNA-Seq data also have the advantage of identifying transcribed regions containing the SNPs, which can then be used to predict loci under selection. Additionally, the identification of SNPs under selection in transcript sequences, coupled with open reading frame (ORF) prediction, can be used to assess how particular SNPs impact the amino-acid replacement rate, which is often associated with adaptation (Ingvarsson 2010; Puillandre et al. 2010; Li et al. 2017; Yu and Li 2018).

The complex life histories of many marine invertebrates and the ocean environment's spatiotemporal variability present significant challenges for estimating biogeographic boundaries based on genetic variation measures (Riginos et al. 2016; Liggins et al. 2019). The rocky intertidal shore is one of the most heterogeneous environments known due to the high seasonality of wave action, as well as rapid changes in salinity and temperature occurring over short periods ranging from months to a few hours (Schaal et al. 2011; Shi et al. 2016; Lathlean et al. 2016). Due to the large effective population sizes found in marine invertebrates, this group may be particularly well-suited to respond to environmental changes (Gagnaire et al. 2015) and is often used as an ecological indicator of environmental disturbances in the oceans (Lima et al. 2011; Sousa et al. 2018).

Littoraria flava (King 1832), a gastropod widely distributed along the Caribbean and Brazilian rock shores, possesses planktotrophic larvae capable of remaining in the water column for 3 to 6 weeks (Reid 1999). Due to their passive dispersal ability, the larvae can cover vast distances, enabling a high level of admixture even among distant populations along the Brazilian coastline (Andrade et al. 2003; Cortez et al. 2021). Although L. flava has high larval dispersal rates, the genetic structure was reported on a microgeographical scale using allozymes and SNPs markers, suggesting the presence of subpopulations within the same rocky shore (Andrade and Solferini 2007; Cortez et al. 2021). In line with these findings, L. flava specimens showed differentially expressed genes in specific sites within the same rocky shore, potentially indicating evidence of local adaptation (Santos et al. 2021).

Based on previous studies indicating signs of local adaptation using differentially expressed genes (Santos et al. 2021), we expect to find a significant correlation between candidate SNPs and marine environmental predictors, suggesting the potential existence of unknown biogeographic barriers and selective forces acting at the molecular level. For that, the present study used genomic resources obtained through GBS and RNA-Seq approaches to identify putative adaptive SNPs in populations from the South West Atlantic Ocean. Using these two different datasets, we were able to explore how selective pressures shape genetic diversity in transcribed and presumed untranscribed regions. This integrative approach allowed us to identify signatures in genomic regions with distinct characteristics. After assessing genetic structure and diversity, potential associations between SNPs and marine environmental predictors were evaluated. Additionally, the signatures of potential selective pressures were investigated by predicting the impact of these SNPs on non-synonymous and noncoding substitutions. In doing so, we were able to explore which potential molecular mechanisms are underlying adaptive responses in the species.

Results

Datasets

After sampling L. flava specimens from 11 locations across the Brazilian coastline (Fig. 1 and supplementary table S1, Supplementary Material online), we obtained 93 specimens for GBS and 40 for RNA-Seq. The library preparation generated 322,479,123 GBS reads. These datasets were originally generated in prior studies to assess both the analysis of differential gene expression through RNA-Seq (Santos et al. 2021) and population structure with GBS (Cortez et al. 2021). After filtering SNPs with minimum allele frequency (MAF) lower than 1% and missing genotypes higher than 20%, 6,298 SNPs from 85 individuals remained (supplementary fig. S1, Supplementary Material online, supplementary table S2, Supplementary Material online). From the RNA-Seq raw data, following methods described in Santos et al. (2021), we obtained 46,863,335 SNPs of 40 different samples. The filtering based on quality (phred score < 30), minimum length of 65 bp, and non-missing data, generated 16,137 SNPs, used in all subsequent analyses 46,863,335 SNPs of 40 different samples (supplementary fig. S1, Supplementary Material online, supplementary table S2, Supplementary Material online). The data used in this study are openly accessible through the SRA database (BioProjects PRJNA665072 and PRJNA656564).

Fig. 1. Littoraria flava sampling details showing both transects and randomly based collection of animals along the Brazilian coastline. The circles indicate the populations where the individuals were sequenced only through GBS technique, the triangles indicate the populations sequenced by both GBS and RNA-Seq and, therefore, sampled along transects. Abbreviations are indicated in supplementary table S1, Supplementary Material online.

An overview of all datasets used in this study, along with the analyses conducted, is presented in Table 1. Further details for each step are provided below.

Table 1 Overview of datasets and corresponding analyses performed in this study

	GBS	RNA-Seq	
	SNPs	Analyses	SNPs	Analyses	
Total	6,298	Selection tests (RDA, LFMM, BayeScan)	16,137	Selection tests (RDA, LFMM, BayeScan)	
Under selection	2,277	…	506	Multivariate analysis (DAPC, PCA); molecular analysis (kN/kS)	
Putative adaptive SNPs	69	Genetic diversity and population structure (θπ, Ho, He, AMOVA, FST); multivariate analysis (DAPC, PCA); functional annotation (BLAST)	4	Functional annotation (BLAST)	
For each GBS or RNA-Seq derived data, the number of SNPs is indicated according to the analysis.

Total, number of total SNPs used; under selection, number of SNPs detected as under selection by at least one selection test (RDA, LFMM or BayeScan) applied; putative adaptive SNPs, number of SNPs commonly identified by at least two out of three selection tests, also represented in Fig. 3.

Putative Candidate SNPs

For both GBS and RNA-Seq, we used two different approaches to identify putative adaptive SNPs: genomic scan based on allele frequencies (BayeScan) and association tests using environmental predictors (latent factor mixed models [LFMMs] and redundancy analysis [RDA]). For the association tests, we first computed the collinearity between each pair of predictors. After removing those with collinearity > 0.8, ten predictors remained: longitude [Long], latitude [Lat], Annual Precipitation [AnnPrec], Precipitation of Wettest Month [PrecWett], Temperature Annual Range [TempAnnRang], Chlorophyll [ChloMean], pH, Salinity, Sea Surface Temperature Range [Sstrange], and Current Velocity Range [Curvssrange] (supplementary fig. S2, Supplementary Material online and supplementary table S3, Supplementary Material online [GBS] and supplementary table S4, Supplementary Material online [RNA-Seq]).

For those 6,298 GBS-derived SNPs, the BayeScan showed that 78 SNPs within 60 loci are potentially under selection (false discovery rate [FDR] < 0.05) (supplementary fig. S3, Supplementary Material online). The RDA full model based on 6,298 SNPs and all ten predictors was significant (P < 0.01) (input as supplementary Appendix S1, Supplementary Material online) showed a correlation of 15% between the environmental and genomic variation. According to the relative arrangement of the samples in the ordination space of axis 1 (P = 0.01), the individual genotypes from Espírito Santo samples (ACF and GAF) might be positively related to the current velocity range (Curvssrange), meaning that individuals from Espírito Santo are under higher current velocity ranges over the year (Fig. 2a). However, these populations present the highest missing data values and, therefore, such outliers might result from differences in genotypes representation.

Fig. 2. Redundancy analysis (RDA) showing the relative contributions of environmental predictors to the genetic structure of putative adaptive genotypes in Littoraria flava. Individual genotypes associations with each predictor are represented for both a) GBS and c) transcriptome candidate SNPs, where SNP are in gray and individuals from distinct locations are represented by different colors according to the map in the left panel. The SNPs significantly associated with each environmental variable are color coded for b) genomic and d) transcriptome dataset (P <0.01). For GBS (a and b), individuals from STF, ACF and GAF are separated from the remaining on the first PC; while for transcriptome (c and d), individuals are geographically separated.

Two candidate SNPs from different loci were associated with Curvssrange (Fig. 2b). The LFMM revealed that 2,266 SNPs were associated with at least one environmental predictor (Fig. 3a). The predictors with the most correlated SNPs were Curvssrange (1,009) followed by pH (830), with 406 SNPs shared between these two. The mapping of all GBS loci containing candidate SNPs resulted in 46 mapped loci (∼5% of 6,298), where about 60% showed the first hits with the gastropod Pomacea canaliculata proteins, which correspond to a pleiotropic regulator, ubiquitin-protein ligase, calcineurin-binding protein, and zinc finger protein (supplementary table S5, Supplementary Material online).

Fig. 3. Venn diagram of the intersection of SNPs containing candidate SNPs identified by RDA, LFMM, and BayeScan for a) GBS dataset of Littoraria flava and b) RNA-Seq. In bold, the final subset of candidate SNPs used for the subsequent analyses.

For the RNA-derived dataset, the BayeScan algorithm identified one single candidate SNP (supplementary fig. S4, Supplementary Material online and Fig. 3b) without any annotation available. The full RDA model based on 16,137 SNPs (input as supplementary Appendix S2, Supplementary Material online) was also significant (P < 0.01) and showed a correlation of 0.2% between the environmental predictors and the genomic variation. A slightly positive relationship of samples from the São João and Anchieta with higher latitude was detected in the first axis (P < 0.01) (Fig. 2c). The sample distribution showed complete segregation of individual genotypes from different locations across the first axis. Also, 135 SNPs within 128 loci were associated with environmental predictors: eight with pH, 31 with Sstrange, 59 with ChloMean, one with Lat, 30 with TempAnnRang, and six with Long (Fig. 2d). None shared two or more predictors. The LFMM analysis detected 374 SNPs associated with at least one environmental predictor, with AnnPrec and Long exhibiting the most associated variants (159 and 150, respectively). Additionally, pH and Curvssrange were associated with 88 and 114 SNPs, respectively. Overall, only four variants were commonly identified by RDA and LFMM (Fig. 3b), which seem related to tRNA ligases and apoptosis (supplementary table S6, Supplementary Material online).

After the adaptive SNPs identification, we looked into functional aspects of those variants. For that, we set the threshold that only variants identified in two of the three analyses would be investigated in terms of potential biological functions. Therefore, the final number was 69 and four candidate SNPs for GBS and transcriptome, respectively (Fig. 3).

Population Diversity and Structure

The measures of genetic diversity and allele polymorphism of the 69 GBS candidate SNPs are presented in Table 2. The θπ values ranged from 0.77 to 7.18 in Gamboa (GAF) and Praia Dura (DF), respectively. The difference between the expected and observed heterozygosity means was significant (P < 0.05), along with a significant and relatively high FIS index (average = 0.606, P < 0.05). When considering all individuals within the same region (Northeast [NE], Southeast [SE], or South [S], see Table 2), the NE exhibited higher diversity despite smaller sampling size. On the other hand, the SE presented the lowest nucleotide diversity, even with the largest sample size. The AMOVA results revealed that most genetic variability (∼63%) occurs within locations (supplementary table S7, Supplementary Material online). The global FST (P = 0.892) and all the pairwise FST values (P > 0.05) failed to reach a level of significance. The discriminant analysis of principal components (DAPC) analysis based on five principal components (PCs), as suggested by the optim.a.score, revealed 11 genetic clusters, which did not present any clear geographic pattern, even when colored according to the regions of Brazilian coastline (Fig. 4 and supplementary fig. S5, Supplementary Material online).

Fig. 4. Discriminant analysis of principal components (DAPC) performed based on 69 GBS candidate SNPs from 85 individuals of Littoraria flava. The distribution of individuals across the two first principal components (PCs) colored according to a) the sampling locations and b) the genetic clusters (K = 11) according to the c) lowest BIC from DAPC, showing that the specimens’ distribution lacks any geographic pattern. Each dot represents an individual. Both DAPCs were performed using five PCs, as suggested by the d) α-score optimization.

Table 2 Littoraria flava genetic diversity indexes based on the 69 GBS candidate SNPs

Location	Samples (n)	θπ	HO	HE	F IS	
Northeast (NE)						
 Region	17	1.7	0.082**	0.341	0.319	
  SBF	11	2.47	0.106**	0.412	0.166	
  ALF	6	5.66	0.120**	0.405	0.615	
Southeast (SE)						
 Region	53	0.628	0.147**	0.339	0.538*	
  ACF	8	1.45	0.125	0.362	1.00**	
  GAF	8	0.77	0.125	0.258	0.00	
  SJF	10	5	0.078**	0.357	1.00*	
  PGF	10	3.27	0.150**	0.327	0.640	
  DF	9	7.18	0.163**	0.422	−0.23	
  ARF	8	3.61	0.125*	0.301	0.740	
South (S)						
 Region	15	0.646	0.066	0.215	1.00*	
  STF	6	2.21	0.055*	0.369	1.00	
  RBF	7	1.16	0.143	0.291	1.00*	
  PIF	2	5	0.556	0.556	NA	
Samples (n)—sampling size; θπ—nucleotide diversity; HO and HE—observed and expected heterozygosity, respectively. The values are shown for each location and for the regions, i.e. all populations from the same region as one single population unit. Abbreviations as in Fig. 1 and supplementary table S1, Supplementary Material online.

*P < 0.05.

**P < 0.01.

Synonymous and Non-Synonymous SNPs

To assess the selection forces acting on the translated amino-acid sequences, we compared the amount of non-synonymous (kN), synonymous (kS) SNPs, and outside coding sequence (oCDS) within neutral and putative adaptive SNPs from the transcriptome dataset only. From all of the 506 adaptive SNPs derived from the transcriptome (according to at least one of the selection analyses, Fig. 3), 303 were annotated within CDS. Of these, 266 were synonymous, and 37 were non-synonymous, resulting in a kN/kS = 0.139 (supplementary table S8, Supplementary Material online). The kN/kS ratio in the 9,287 neutral SNPs annotated within CDS was 0.148 (supplementary table S8, Supplementary Material online). Pearson's Chi-square test showed no difference in the kN/kS proportion between the set of putatively adaptive SNPs and the putatively neutral ones (P > 0.7). Also, none of the four SNPs commonly found by RDA and LFFM (putative candidate SNPs section) were non-synonymous changes (Fig. 5).

Fig. 5. Contribution of the top 200 SNPs in the descending contribution rank to the discriminant functions 1 a) and 2 b) of discriminant analysis of principal components (DAPC) and principal components 1 c) and 2 d) of principal components analysis (PCA). Each bar corresponds to the percentage of contribution of the SNP to the component analysis axis. The yellow bars represent the SNPs identified by RDA and LFMM, all of which were synonymous SNPs. At the first iteration, the top ten SNPs from the descending rank of were tested and, with each iteration, the following SNP of the rank is added to the test. The dashed red line represents the P-value log (0.05). Red asterisks are P < 0.05.

Subsequently, we described the differences in adaptive SNP frequency between localities by applying both principal component analysis (PCA) and DAPC from R packages “FactorMineR” and “adegenet”. We also evaluated the contribution of kN SNPs to the obtained principal components. These analyses revealed that it was possible to ascertain the role of substitutions affecting the translated amino-acid sequence on local selection. The first and second axis of the PCA represented 56.9% and 27.1% of the variance, respectively (supplementary fig. S6, Supplementary Material online), totaling 84% of all variation within the data. The DAPC was performed with six PCs as suggested by the function optimal.a.score (supplementary fig. S6, Supplementary Material online). The eigenvalues of the first and second discriminant functions, interpreted as the variance ratio between the defined cluster to the variance within these clusters, were 181.2 and 8.4, indicating a crucial role of the first component in segregating the localities. Both the PCA and DAPC resulted in the first component separating AR and DF from SJF and ACF. Moreover, the second axis separates ACF from SJF (supplementary fig. S6, Supplementary Material online). Neither the first PCA component nor the first DAPC discriminant function displayed a significant excess of kN SNPs at the top of the contribution ranking as confirmed by the Wilcoxon rank sum test (Fig. 5). Conversely, a nonrandom (P < 0.05) distribution of kN, kS, and oCDS SNPs in the contribution rank was only observed when testing the top 94 and 98 SNPs from the second principal component and second discriminant function (Fig. 5).

Discussion

We combined genome and transcriptome-wide SNPs obtained from 125 specimens with seascape attributes to detect potential local adaptation in a highly dispersive marine invertebrate. While previous neutral genetic structures indicated the presence of three genetic clusters (Cortez et al. 2021), we found 11 putative adaptive genetic clusters without any clear geographic pattern. In addition, we observed different segregation patterns across localities for adaptive SNPs derived from GBS and RNA-Seq, revealing that transcribed and supposedly non-transcribed regions seem to be experiencing distinct selective forces. Remarkably, non-synonymous substitutions did not contribute significantly within the transcribed regions, in contrast to the highly relevant role of synonymous and oCDS substitutions for genetic structure.

Here, the implemented methods showed good performance for both datasets. Despite the difference in the number of identified SNPs, we were able to recover SNPs commonly found in at least two out of the three analyses, reducing the false positives and providing a more robust dataset. For the association tests, even though public databases do not yet encompass all conditions capable of predicting genetic diversity, the currently available seascape attributes have provided relevant information about biogeographic barriers and selective forces in the marine realm (Riginos et al. 2016), particularly when integrated with high-throughput sequencing technologies (Galindo et al. 2010; Westram et al. 2014). We have found divergent results with RNA and GBS-derived SNPs, suggesting different independent selective pressures acting in transcribed and supposedly non-transcribed genomic regions (Pespeni et al. 2012). Additionally, only transcriptome candidate SNPs showed a clear geographic distribution pattern of genotypes in the RDA (Fig. 2c), probably reflecting environmental differences, such as human activities, in these areas (Sterza and Fernandes 2006). Nonetheless, the distribution of individual genotypes in RDA ordination can also be influenced by missing data rate. For example, populations ACF and GAF show higher rates of missing data in the GBS dataset, and their specimens are depicted sparsed in the RDA plot (Fig. 2). The pattern observed could be related to the missing data, so these populations’ results are interpreted carefully.

Eleven adaptive genetic clusters were found when analyzing putative adaptive SNPs from GBS, despite lacking any geographic pattern (Fig. 4). Previous findings with putatively neutral SNPs showed three discrete genetic clusters from the same collection locations, presented in Cortez et al. (2021). Even presenting a high level of admixture, the marine currents seem to be able to shape the distribution of genomic variation of L. flava and other invertebrates (Benestan et al. 2016; Miller et al. 2019). The candidate SNPs still showed a nonsignificant FST and high nucleotide diversity. Also, the diversity indices were higher for populations in Northeastern and Southern localities, despite the smallest sampling size when considering both locations or regions as a population unit (Table 2). For the Southeastern localities, besides the excess of heterozygotes for some local populations, a higher admixture level or a common origin of the individual genotypes is supported. Most genomic variation is mainly explained by individuals within the rocky shores instead of in large-scale coastal regions. Together with the AMOVA and FST results from neutral SNPs (Cortez et al. 2021), these findings keep supporting highly connected L. flava populations across the Brazilian rocky shore, even though unclear sources of selective forces might affect the genetic diversity distribution. Thus, although the current velocity range was the only predictor commonly found for the GBS and transcriptome candidate SNPs, other predictors, including pH, water surface temperature, chlorophyll, precipitation and temperature range, influence the distribution of transcriptome-derived SNPs. These conditions are known to be relevant for the performance and survival of species that inhabit rocky shores and are exposed to both marine and terrestrial conditions (Harley and Helmuth 2003; Fuchs et al. 2010; Benestan et al. 2016; Miller et al. 2019; Sauerland et al. 2019; Takeuchi et al. 2020; Mendes et al. 2022).

From RNA-Seq candidate SNPs, two loci gene products are related to apoptosis: serine/threonine-protein phosphatase 6 (PP6) and apoptosis-inducing factor 1 (AIF1) (supplementary table S6, Supplementary Material online). AIF1 was also identified and characterized in Crassostrea gigas by Qiao et al. (2021) as participating in apoptosis triggered by immune responses. Apoptosis may be triggered during immune responses caused by environmental disturbances, such as pollution, hormones, temperature, and pH variations (Romero et al. 2015; Qiao et al. 2021). In gastropods, cell death has also been reported to follow cellular damage after exposure to heavy metals and xenobiotics (Chabicovsky et al. 2004; Russo and Madec 2007). A PP6 loci variant was identified as under selection, which was also found up-regulated in São João gastropods (Santos et al. 2021), an area known for fishing, boat activities, and high levels of pollutants in the water, which may lead to a stressful scenario for fauna living along rocky shores (Oigman-Pszczol and Creed 2007; do Nascimento et al. 2018). Interestingly, PP6 is believed to be involved in shell formation and biomineralization processes and ocean acidification changes can represent a higher vulnerability for gastropods larvae development or calcification (Kurihara 2008; Ellis et al. 2009; Gazeau et al. 2010; Duarte et al. 2013; Rajan and Vengatesen 2020).

We investigated the impact of local selection on gene variation. For this purpose, transcriptomic data allowed us to compare the contribution of kN SNPs, kS SNPs, and oCDS SNPs to potential local adaptation. We found that the ratio of kN/kS was similar between putatively adaptive and the whole set of SNPs (supplementary table S8, Supplementary Material online). This similarity suggests these variants may play equivalent roles or have no adaptive significance. Accordingly, the Wilcoxon rank sum tests for the first component of both the PCA and DAPC, which segregates ARF and PGF from SJF and ACF (supplementary fig. S6, Supplementary Material online), suggests that genetic clustering between these localities did not arise mainly from the natural selection of polymorphic proteins, but from significant contributions of synonymous SNPs and oCDS SNPs. On the other hand, Wilcoxon rank sum tests indicates that kN SNPs are not randomly distributed through the contribution rank for the second axis of PCA and DAPC (P < 0.05) (Fig. 5), both of which segregate SJF from ACF (supplementary fig. S6, Supplementary Material online). The eight kN SNPs present in the significant Wilcoxon tests have a high contribution to the segregation of the aforementioned localities, which are close to each other. The nonrandom ranking of kN SNPs contribution is suggestive of amino-acid level divergences between localities, but there is an enormous contribution of kS SNPs and oCDS SNPs for all components. Further nonrandom distribution of kN SNPs in the contribution rank can hardly be confirmed by visual inspection of the data (Fig. 5).

Through mechanisms previously suggested, such as degradation rates and splicing efficiency of the mRNA (Chamary et al. 2006; Liu et al. 2021), these SNPs might be affecting the expression patterns between localities as a result of local selection. Synonymous SNPs under selection affecting translation efficiency were detected in human evolution (Waldman et al. 2011). An alternative hypothesis is that these SNPs are not directly selected but are genetically linked to polymorphic promoters, usually close to their respective transcript in eukaryotic genomes (Haberle and Stark 2018). Both hypotheses could participate in the differentiation of gene expression across localities, a pattern observed for L. flava (Santos et al. 2021). However, linkage between kS and kN SNPs present in unsampled transcripts should not be discarded. (P < 0.05)

Our findings suggest that amino-acid polymorphisms (resulting from kN SNPs) play a minor role in local selection for L. flava adults. SNPs selected because of different ecological conditions were mostly kS and oCDS. Presence of putatively adaptive kS SNPs was previously observed in oysters from the genus Crassostrea (Zhong et al. 2014) and in fish from the species Merluccius merluccius (Milano et al. 2014), and positively selected sites in UTRs were found in Drosophila species (Haddrill et al. 2008). The contribution of adaptive oCDS SNPs to genetic structure was discussed for populations of sea urchins, although these were in regions upstream to the transcription start site, not in UTRs, as in our results (Pespeni et al. 2012). In line with that, reduction of spines in threespine sticklebacks provides an example of adaptive genetic variation in untranscribed regions (Shapiro et al. 2004; Chan et al. 2010). Despite the presented evidence for kS and oCDS SNPs participation in local selection, some authors treat kS SNPs as false positives instead of considering the possibility of selection through other means (Cullingham et al. 2014). Combined with previous researches, our results suggest that these SNPs should not be overlooked as they represent a substantial portion of the polymorphisms putatively under selection.

Our study presents an attempt to understand the role of kS and transcripted oCDS SNPs in local selection in a non-model organism. These findings have great potential for future local adaptation studies regarding the mechanisms driving adaptive evolution in natural populations. Furthermore, implementing molecular research alongside environmental association tests highlights the benefits of integrative approaches for understanding natural selection dynamics, especially for the complex marine realm.

Materials and Methods

Biological Sampling

We collected L. flava specimens from rocky shores at 11 locations along the Brazilian coastline (supplementary table S1, Supplementary Material online) under the Instituto Chico Mendes de Conservação da Biodiversidade license no. 56726-1. In four of these localities (Araçá, Praia Dura, São João, and Anchieta), samples were collected along horizontal transects to examine the effect of microhabitat distribution on SNP frequencies and gene expression (Fig. 1). The transect on the rocky shore extended from the furthest point from the sea to the closest point to the sea where the species was present, varying from 0 to 64 m of distance. Sampling methodology details can be found in Cortez et al. (2021) and Santos et al. (2021). Whole-body individuals were then separately subjected to RNA and DNA extraction.

GBS and SNPs Call

We extracted genomic DNA following Doyle and Doyle's (1987) protocol and quantified DNA concentrations using a dsDNA BR Assay kit (Invitrogen, Waltham, MA, USA) with a Qubit v3 fluorometer (ThermoFisher Scientific, Waltham, MA, USA). Libraries were constructed based on the GBS method following the protocol of Elshire et al. (2011) and using the PstI restriction enzyme (5′-CTGCAG-3′). Once ligated to the barcode and common adaptors, the products were PCR amplified using generic primers matching the common adaptors under the following conditions: 5 min at 72 °C, 30 s at 98 °C, 18 cycles of 10 s at 98 °C, 30 s at 65 °C, and 30 s at 72 °C, and an extension step of 5 min at 72 °C. The size and quality of the DNA fragments were confirmed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) with the Agilent DNA 1000 kit and by qPCR on a Light Cycler 480II (Roche). The Kapa Biosystems kit was used for library quantitation. The GBS libraries were constructed at the EcoMol Consultoria (Piracicaba, SP, Brazil) and sequenced at the Center for Functional Genomics Applied to Agriculture and Agroenergy (Animal Biotechnology Laboratory, LZT/ESALQ/USP, Piracicaba, SP, Brazil) on a HiSeq 2500 platform (Illumina Inc., San Diego, CA, USA).

The quality of GBS raw reads was checked with FastQC v.0.11.8 (Andrews 2010). Seqyclean v.1.9.9 (Zhbannikov et al. 2017) through the UniVec database (NCBI, ftp://ftp.ncbi.nlm.nih.gov/pub/UniVec/) was utilized to remove adapter sequences, vectors, oligonucleotides, and reads with average Phred (QS) quality below 24 and smaller than 50 bp. Individual read assignments, paralog identification, and read clustering into consensus sequence for each locus were performed with the iPyrad v.0.7.28 program (Eaton 2014). All program-specific input formats were obtained with PGDSpider v.2.1.15 (Lischer and Excoffier 2012). PLINK v.2.0 (Purcell et al. 2007) was used to remove SNPs with a MAF of <0.01 and missing genotypes per SNP of >0.20. Details on the filtering results can be found in Cortez et al. (2021).

RNA-Seq, Transcriptome Assembly, and SNPs Call

We isolated the RNA from entire L. flava animals obtained at each site using Trizol, according to Riesgo et al. (2012). The RNA quantity was verified using the Qubit fluorometer and NanoDrop spectrophotometer (ThermoFisher Scientific). The integrity of the samples was confirmed using BioAnalyzer (Agilent Technologies Inc.), and only samples with RNA integrity number (RIN) values of >6.0 were selected for subsequent analysis. A total of 40 libraries were constructed and sequenced as described in Santos et al. (2021).

After sequencing, we visualized the quality of the raw data generated with FastQC software. All the reads for adapters, primers, and quality were filtered using SeqyClean with QS with the mean and edge minimum score values set at 24 and 30, respectively, and a minimum length of 65 bp. We performed a de novo assembly using Trinity v.2.8.4 (Grabherr et al. 2011), as described in Santos et al. (2021) and used the TransDecoder package (http://transdecoder.sourceforge.net/) to identify the candidate coding regions. Next, the functional annotation was performed using the Trinotate pipeline (https://trinotate.github.io/) with BLASTx in the following databases: UniProt (uniref90 + SwissProt) with a cutoff value of 1e10-3; Gene Ontology (GO) (Ashburner et al. 2000), with GO terms: biological process (BP), molecular function (MF), and cellular component (CC); and KEGG (Kyoto Encyclopedia of Genes and Genomes) (Kanehisa et al. 2012).

We then mapped the reads of the 40 L. flava samples against the reference unigenes with BLASTx hits using Bowtie2 v.2.3.4.3 (Langmead and Salzberg 2012). Using the mpileup command of Samtools v.1.3 (Li et al. 2009), we detected SNPs on annotated transcripts using the parameters -q20, -C50, -A (do not discard anomalous read pairs) and -B (disable base alignment quality inference). With BCFtools v.1.3 (Li et al. 2009) for SNP calling, we excluded variant base quality of <30 and the sum of the coverage of reads with alternative alleles in the forward and reverse (DP4) strands of <10 from the analyses to avoid artifacts. We only considered variants present in at least four individuals (10%), with a minimum frequency of 98% and MAF < 0.01 for downstream analysis, using the options –mac 4, --max-missing 0.98 and --maf 0.01in vcftools v.0.1.16 (Danecek et al. 2011).

Putative Candidate SNPs

We conducted a genome scan and two environmental association tests to identify the best candidate SNPs for local adaptation using RNA-Seq and GBS read datasets. By employing the Bayesian approach from BayeScan v.2.0 (Foll and Gaggiotti 2008), we used a genome scan to detect deviations from neutrality. We set a prior odd of 10, which assumes that the neutral model is ten times more likely than the selection model. The program ran 50,000 iterations followed by 500,000 simulations, with a burn-in of 1% and a FDR of 5%.

We performed environmental association tests using RDA (Forester et al. 2018) and LFMM (Frichot and François 2015). Previous studies have suggested that several environmental variables, such as salinity (She et al. 2018), sea temperature (Dwane et al. 2023), and dissolved oxygen (San et al. 2021), can drive local adaptation in marine invertebrate species. Therefore, we downloaded continental and oceanic environmental data from sampling locations from WorldClim (www.worldclim.org) and Bio-ORACLE databases (Tyberghein et al. 2012). We extracted data for each locality using the sdmpredictors library (Bosch et al. 2017) from the statistical software R v.3.6 (R Core Team 2013). We tested the collinearity of data using RDA from the package “vegan” 2.6.2 (Oksanen et al. 2013). We removed all values greater than 0.8, retaining only ten predictors (Long, Lat, AnnPrec, PrecWett, TempAnnRang, ChloMean, pH, Salinity, Sstrange, Curvssrange; see Results section). We tested the significance of the full model and each constrained axis using an analysis of variance with 1,000 permutations. We applied a standard deviation cutoff of ±3 from the mean SNP loading for each axis (Forester et al. 2018).

We implemented the LFMM test using R package “LEA” 3.15 (Frichot and François 2015) for the same ten environmental predictors used in RDA. We tested each predictor for each K using a burn-in of 2,000 followed by 20,000 iterations, setting the number of latent factors according to the minimal cross-entropy criterion of the snmf function (K varying from 1 to 20). We adjusted the P-values based on the median z-score (Storfer et al. 2018) and adopted a significance level of 1% (P < 0.01). We defined the set of candidate SNPs for local adaptation as those found in at least two of the three analyses (i.e. BayeScan, RDA, and LFMM) for both datasets, which are hereafter referred to as GBS and transcriptome candidate SNPs.

The GBS candidate SNPs were mapped against the previously annotated transcriptome of L. flava using Bowtie2 v.2.4.4 (Langmead and Salzberg 2012). Functional annotation for transcriptome candidate markers was performed with the BLASTn Search tool from BLAST v.2.9.0 (Camacho et al. 2009) using the nucleotide collection and e-value threshold of 10−3.

Genetic Diversity and Structure for Candidate SNPs

For the GBS candidate SNPs, we calculated multi-locus expected heterozygosity (HE), observed heterozygosity (HO), nucleotide diversity (θπ), and the fixation index FIS (Weir and Cockerham 1984) for the GBS candidate SNPs using Arlequin v.3.5 (Excoffier and Lischer 2010). The Bartlett test was applied for the variance of observed and expected heterozygosity and obtained the significance of the mean differences using the paired t-test. We tested the significance level with 10,000 permutations (P < 0.05). The genetic differentiation was evaluated using the AMOVA from Arlequin, considering (i) each locality and (ii) each region (NE, SE, and S, see supplementary table S1, Supplementary Material online), based on the unbiased FST estimator θ (Weir and Cockerham 1984).

We used DAPC (Jombart et al. 2010) as an assumption-free genetic clustering method to assess adaptive genetic structure. This approach was preferred since the genetic clustering algorithm from STRUCTURE (Pritchard et al. 2000) requires Hardy–Weinberg equilibrium and lack of linkage disequilibrium in ancestral populations (Frichot et al. 2014). We performed the DAPC using the optimal number of PCs (n = 5, see results) suggested by the optim.a.score function from the “adegenet” R package (Jombart 2008) and defined genetic clusters as the locality of the individuals using the find.cluster function.

Synonymous and Non-Synonymous SNPs Analyses

To understand how local selection affects SNP frequency in transcribed sequences, we annotated and classified all filtered SNPs detected in the RNA-Seq data as synonymous, non-synonymous, or oCDS using SnpEff v.5.0 (Cingolani et al. 2012) with the candidate CDS previously identified by TransDecoder as a reference.

We filtered out transcripts containing more than one predicted ORF since these could result in ambiguously annotated SNPs. As only 46 GBS SNPs putatively under selection could be mapped to transcribed regions, these SNPs were not considered here. The number of SNPs was counted and classified as non-synonymous (kN) (except for loss of sense), synonymous (kS), and untranslated regions (i.e. oCDS) for both neutral and putatively adaptive SNPs. For that, we used a custom bash script (supplementary Appendix S3, Supplementary Material online) and a variant calling file (vcf) generated by SnpEff. We defined neutral SNPs as those identified by none of the analyses (BayeScan, RDA, and LFMM). As too few SNPs were found by at least two predictors, for this analysis, we considered putatively adaptive SNPs as those 506 SNPs identified by at least one of the three analyses. Pearson's Chi-squared test using the R package “stats” (R Core Team 2013) was performed to assess the similarity in the proportion of kN and kS SNPs between the two sets.

To describe the contribution to local selection provided by SNPs that change the translated amino-acid chain (kN) and those that do not (kS and oCDS), we used principal component approaches. We performed a PCA using the “FactoMineR” (Lê et al. 2008) package and a DAPC using the “adegenet” package. The PCA was calculated using each locality as an observation and bi-allelic amino-acid polymorphism frequencies as variables. We set each locality as a genetic cluster for DAPC (PCs = six, according to the optim.a.score function). The PCA and DAPC components segregated individuals from different localities based on SNPs’ frequency variance. We expect that variation in SNP frequencies across different localities may result from distinct selective forces reflecting the particularities of the biotic and abiotic conditions. These dynamics could, then, lead to local adaptation.

Afterward, we compared the contribution of kN SNPs to the contribution of kS and oCDS SNPs (previously annotated by snpEff) in the first two components of the PCA and DAPC. For that, we ranked the SNPs by their contribution to each component in descending order and submitted the top 200 SNPs (∼90% of total contribution) of each component to a one-tailed Wilcoxon rank sum test (Bauer 1972). We compared the ranks of kN SNPs to the ranks of kS and oCDS SNPs. With this test, we assessed if kN SNPs are randomly distributed through the contribution rank or if they are more prone to be at the top. We used a custom R script (supplementary Appendix S4, Supplementary Material online) that starts with the ten top-ranked SNPs and successively adds the following SNPs in the descending contribution, applying a one-tailed Wilcoxon rank sum test at each iteration, until the contribution of the top 200 SNPs was tested. The iterations increase the sensibility by removing SNPs with a negligible contribution to the components, which could add noise to the analyses.

Supplementary Material

evae194_Supplementary_Data

Acknowledgments

We are thankful to the Centro de Biologia Marinha da Universidade de São Paulo (CEBIMar) staff, the EcoMol staff for the GBS library preparation, the Centro de Genômica Funcional do Laboratório de Biotecnologia Animal (ESALQ-USP) for the RNA-Seq sequencing, and the Darwin server administrators for all the assistance. This work was supported by the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP, processes 2015/20139-9 and 2018/05118-3) and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq, process 155636/2018-9).

Supplementary Material

Supplementary material is available at Genome Biology and Evolution online.

Author Contributions

T.C., C.A.S., and S.C.S.A. designed the study. T.C. conducted all selection tests and GBS-derived and population analyses. C.A.S. performed all RNA-Seq bioinformatics, annotation, and mapping analyses. G.G.S. performed molecular evolutionary analyses. All contributed to the writing and editing of the manuscript. All authors read and approved the final manuscript.

Conflict of Interest

The authors have no conflicts of interest to declare.

Data Availability

RNA-Seq and GBS-derived raw sequence reads are deposited in the SRA (BioProjects PRJNA665072 and PRJNA656564, respectively).
==== Refs
Literature Cited

Andrade  SCS, Magalhães  CA, Solferini  VN. Patterns of genetic variability in Brazilian littorinids (Mollusca): a macrogeographic approach. J Zool Syst Evol Res. 2003:41 (4 ):249–255. 10.1046/j.1439-0469.2003.00227.x.
Andrade  SCS, Solferini  VN. Fine-scale genetic structure overrides macro-scale structure in a marine snail: nonrandom recruitment, demographic events or selection?  Biol J Linnean Soc. 2007:91 (1 ):23–36. 10.1111/j.1095-8312.2007.00782.x.
Andrews  S.  2010. FastQC: a quality control tool for high throughput sequence data. Accessed 2023 Oct 10. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
Ashburner  M, Ball  CA, Blake  JA, Botstein  D, Butler  H, Cherry  JM, Davis  AP, Dolinski  K, Dwight  SS, Eppig  JT ,  et al  Gene ontology: tool for the unification of biology. The gene ontology consortium. Nat Genet. 2000:25 (1 ):25–29. 10.1038/75556.10802651
Bauer  DF . Constructing confidence sets using rank statistics. J Am Stat Assoc. 1972:67 (339 ):687–690. 10.1080/01621459.1972.10481279.
Bay  RA, Palumbi  SR. Multilocus adaptation associated with heat resistance in reef-building corals. Curr Biol. 2014:24 (24 ):2952–2956. 10.1016/j.cub.2014.10.044.25454780
Benestan  L, Quinn  BK, Maaroufi  H, Laporte  M, Clark  FK, Greenwood  SJ, Rochette  R, Bernatchez  L. Seascape genomics provides evidence for thermal adaptation and current-mediated population structure in American lobster (Homarus americanus). Mol Ecol. 2016:25 (20 ):5073–5092. 10.1111/mec.13811.27543860
Bosch  S, Tyberghein  L, De Clerck  O. sdmpredictors: An R package for species distribution modelling predictor datasets. Marine Species Distributions: From data to predictive models, 2017.
Camacho  C, Coulouris  G, Avagyan  V, Ma  N, Papadopoulos  J, Bealer  K, Madden  TL. BLAST+: architecture and applications. BMC Bioinformatics. 2009:10 (1 ):421. 10.1186/1471-2105-10-421.20003500
Chabicovsky  M, Klepal  W, Dallinger  R. Mechanisms of cadmium toxicity in terrestrial pulmonates: programmed cell death and metallothionein overload. Environ Toxicol Chem. 2004:23 (3 ):648–655. 10.1897/02-617.15285358
Chamary  JV, Parmley  JL, Hurst  LD. Hearing silence: non-neutral evolution at synonymous sites in mammals. Nat Reviews Gen. 2006:7 (2 ):98–108. 10.1038/nrg1770.
Chan  YF, Marks  ME, Jones  FC, Villarreal  G  Jr, Shapiro  MD, Brady  SD, Southwick  AM, Absher  DM, Grimwood  J, Schmutz  J ,  et al  Adaptive evolution of pelvic reduction in sticklebacks by recurrent deletion of a Pitx1 enhancer. Science. 2010:327 (5963 ):302–305. 10.1126/science.1182213.20007865
Cingolani  P, Platts  A, Wang  LL, Coon  M, Nguyen  T, Wang  L, Land  SJ, Lu  X, Ruden  DM. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin).  2012:6 (2 ):80–92. 10.4161/fly.19695.22728672
Cortez  T, Amaral  RV, Sobral-Souza  T, Andrade  SC. Genome-wide assessment elucidates connectivity and the evolutionary history of the highly dispersive marine invertebrate Littoraria flava (Littorinidae: Gastropoda). Biol J Linnean Soc. 2021:133 (4 ):999–1015. 10.1093/biolinnean/blab055.
Cullingham  CI, Cooke  JE, Coltman  DW. Cross-species outlier detection reveals different evolutionary pressures between sister species. New Phytol. 2014:204 (1 ):215–229. 10.1111/nph.12896.24942459
Danecek  P, Auton  A, Abecasis  G, Albers  CA, Banks  E, DePristo  MA, Handsaker  RE, Lunter  G, Marth  GT, Sherry  ST, et al  The variant call format and VCFtools. Bioinformatics. 2011:27 (15 ):2156–2158. 10.1093/bioinformatics/btr330.21653522
do Nascimento  MTL, de Oliveira Santos  AD, Felix  LC, Gomes  G, de Oliveira E Sá  M, da Cunha  DL, Vieira  N, Hauser-Davis  RA, Baptista Neto  JA, Bila  DM. Determination of water quality, toxicity, and estrogenic activity in a nearshore marine environment in Rio de Janeiro, Southeastern Brazil. Ecotoxicol Environ Saf. 2018:149 :197–202. 10.1016/j.ecoenv.2017.11.045.29175346
Doyle  JJ, Doyle  JL. A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochem Bull. 1987:19 :11–15.
Duarte  C, Hendriks  IE, Moore  TS, Olsen  Y, Steckbauer  A, Ramajo  L, Carstensen  J, Trotter  J, McCulloch  M. Is ocean acidification an open-ocean syndrome? Understanding anthropogenic impacts on seawater pH. Estuaries Coast. 2013:36 (2 ):221–236. 10.1007/s12237-013-9594-3.
Dwane  C, Rezende  EL, Tills  O, Galindo  J, Rolán-Alvarez  E, Rundle  S, Truebano  M. Thermodynamic effects drive countergradient responses in the thermal performance of Littorina saxatilis across latitude. Sci Total Environ. 2023:863 :160877. 10.1016/j.scitotenv.2022.160877.36521622
Eaton  DA . PyRAD: assembly of de novo RADseq loci for phylogenetic analyses. Bioinformatics. 2014:30 (13 ):1844–1849. 10.1093/bioinformatics/btu121.24603985
Ellis  RP, Bersey  J, Rundle  SD, Hall-Spencer  JM, Spicer  JI. Subtle but significant effects of CO2 acidified seawater on embryos of the intertidal snail, Littorina obtusata. Aquat Biol. 2009:5 :41–48. 10.3354/ab00118.
Elshire  RJ, Glaubitz  JC, Sun  Q, Poland  JA, Kawamoto  K, Buckler  ES, Mitchell  SE. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS One. 2011:6 (5 ):e19379. 10.1371/journal.pone.0019379.21573248
Excoffier  L, Lischer  HE. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010:10 (3 ):564–567. 10.1111/j.1755-0998.2010.02847.x.21565059
Foll  M, Gaggiotti  O. A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. Genetics. 2008:180 (2 ):977–993. 10.1534/genetics.108.092221.18780740
Forester  BR, Lasky  JR, Wagner  HH, Urban  DL. Comparing methods for detecting multilocus adaptation with multivariate genotype–environment associations. Mol Ecol. 2018:27 (9 ):2215–2233. 10.1111/mec.14584.29633402
Frichot  E, François  O. LEA: an R package for landscape and ecological association studies. Methods Ecol Evol. 2015:6 (8 ):925–929. 10.1111/2041-210X.12382.
Frichot  E, Mathieu  F, Trouillon  T, Bouchard  G, François  O. Fast and efficient estimation of individual ancestry coefficients. Genetics. 2014:196 (4 ):973–983. 10.1534/genetics.113.160572.24496008
Fuchs  HL, Solow  AR, Mullineaux  LS. Larval responses to turbulence and temperature in a tidal inlet: habitat selection by dispersing gastropods?  J Mar Res. 2010:68 (1 ):153–188. 10.1357/002224010793079013.
Gagnaire  P-A, Broquet  T, Aurelle  D, Viard  F, Souissi  A, Bonhomme  F, Arnaud-Haond  S, Bierne  N. Using neutral, selected, and hitchhiker loci to assess connectivity of marine populations in the genomic era. Evol Appl. 2015:8 (8 ):769–786. 10.1111/eva.12288.26366195
Gagnaire  P-A, Normandeau  E, Bernatchez  L. Comparative genomics reveals adaptive protein evolution and a possible cytonuclear incompatibility between European and American eels. Mol Biol Evol. 2012:29 (10 ):2909–2919. 10.1093/molbev/mss076.22362081
Galindo  J, Grahame  JW, Butlin  RK. An EST-based genome scan using 454 sequencing in the marine snail Littorina saxatilis. J Evol Biol. 2010:23 (9 ):2004–2016. 10.1111/j.1420-9101.2010.02071.x.20695960
Gazeau  F, Gattuso  JP, Dawber  C, Pronker  AE, Peene  F, Peene  J, Heip  CHR, Middelburg  JJ. Effect of ocean acidification on the early life stages of the blue mussel Mytilus edulis. Biogeosciences. 2010:7 (7 ):2051–2060. 10.5194/bg-7-2051-2010.
Gleason  LU, Burton  RS. RNA-Seq reveals regional differences in transcriptome response to heat stress in the marine snail Chlorostoma funebralis. Mol Ecol. 2015:24 (3 ):610–627. 10.1111/mec.13047.25524431
Grabherr  MG, Haas  BJ, Yassour  M, Levin  JZ, Thompson  DA, Amit  I, Adiconis  X, Fan  L, Raychowdhury  R, Zeng  Q ,  et al  Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011:29 (7 ):644–652. 10.1038/nbt.1883.21572440
Haberle  V, Stark  A. Eukaryotic core promoters and the functional basis of transcription initiation. Nat Rev Mol Cell Biol. 2018:19 (10 ):621–637. 10.1038/s41580-018-0028-8.29946135
Haddrill  PR, Bachtrog  D, Andolfatto  P. Positive and negative selection on noncoding DNA in Drosophila simulans. Mol Biol Evol. 2008:25 (9 ):1825–1834. 10.1093/molbev/msn125.18515263
Harley  CD, Helmuth  BS. Local- and regional-scale effects of wave exposure, thermal stress, and absolute versus effective shore level on patterns of intertidal zonation. Limnol Oceanogr. 2003:48 (4 ):1498–1508. 10.4319/lo.2003.48.4.1498.
Huang  Y, Chain  FJ, Panchal  M, Eizaguirre  C, Kalbe  M, Lenz  TL, Samonte  IE, Stoll  M, Bornberg-Bauer  E, Reusch  TBH ,  et al  Transcriptome profiling of immune tissues reveals habitat-specific gene expression between lake and river sticklebacks. Mol Ecol. 2016:25 (4 ):943–958. 10.1111/mec.13520.26749022
Ingvarsson  PK . Natural selection on synonymous and nonsynonymous mutations shapes patterns of polymorphism in Populus tremula. Mol Biol Evol. 2010:27 (3 ):650–660. 10.1093/molbev/msp255.19837657
Jombart  T . adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008:24 (11 ):1403–1405. 10.1093/bioinformatics/btn129.18397895
Jombart  T, Devillard  S, Balloux  F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 2010:11 (1 ):94. 10.1186/1471-2156-11-94.20950446
Kanehisa  M, Goto  S, Sato  Y, Furumichi  M, Tanabe  M. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Res. 2012:40 (D1 ):D109–D114. 10.1093/nar/gkr988.22080510
King  P, Broderip  W. Description of the Cirrhipeda, Conchifera and Mollusca: in a collection formed by the Officers of HMS Adventure and Beagle employed between the years 1826 and 1830 in surveying the southern coasts of South America: including the Straits of Magalhaens and the Coast of Tierra Del Fuego. J Zool. 1832:5 :332–349. 10.1093/nar/gkr988.
Kurihara  H . Effects of CO2-driven ocean acidification on the early developmental stages of invertebrates. Mar Ecol Prog Ser. 2008:373 :275–284. 10.3354/meps07802.
Langmead  B, Salzberg  SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012:9 (4 ):357–359. 10.1038/nmeth.1923.22388286
Lathlean  JA, Seuront  L, McQuaid  CD, Ng  TP, Zardi  GI, Nicastro  KR. Size and position (sometimes) matter: small-scale patterns of heat stress associated with two co-occurring mussels with different thermoregulatory behaviour. Mar Biol. 2016:163 (9 ):1–11. 10.1007/s00227-016-2966-z.
Lê  S, Josse  J, Husson  F. FactoMineR: an R package for multivariate analysis. J Stat Softw. 2008:25 (1 ):1–18. 10.18637/jss.v025.i01.
Li  H, Handsaker  B, Wysoker  A, Fennell  T, Ruan  J, Homer  N, Marth  G, Abecasis  G, Durbin  R; 1000 Genome Project Data Processing Subgroup. The sequence alignment/map format and SAMtools. Bioinformatics. 2009:25 (16 ):2078–2079. 10.1093/bioinformatics/btp352.19505943
Li  X, Shi  L, Zhou  Y, Xie  H, Dai  X, Li  R, Chen  Y, Wang  H. Molecular evolutionary mechanisms driving functional diversification of α-glucosidase in Lepidoptera. Sci Rep. 2017:7 (1 ):45787. 10.1038/srep45787.28401928
Liggins  L, Treml  EA, Riginos  C. Seascape Genomics: Contextualizing Adaptive and Neutral Genomic Variation in the Ocean Environment. In: Oleksiak  M, Rajora  O, editors. Population Genomics: Marine Organisms. Springer; 2019. p. 171–218. 10.1007/13836_2019_68.
Lima  D, Reis-Henriques  MA, Silva  R, Santos  AI, Castro  LFC, Santos  MM. Tributyltin-induced imposex in marine gastropods involves tissue-specific modulation of the retinoid X receptor. Aquat Toxicol. 2011:101 (1 ):221–227. 10.1016/j.aquatox.2010.09.022.21036407
Lischer  HEL, Excoffier  L. PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics. 2012:28 (2 ):298–299. 10.1093/bioinformatics/btr642.22110245
Liu  Y, Yang  Q, Zhao  F. Synonymous but not silent: the codon usage code for gene expression and protein folding. Annu Rev Biochem. 2021:90 (1 ):375–401. 10.1146/annurev-biochem-071320-112701.33441035
Mendes  CB, Cortez  T, Santos  CS, Sobral-Souza  T, Santos  AD, Sasaki  DK, Augusto Silva  D, Dottori  M, Andrade  SCS. Seascape genetics in a polychaete worm: disentangling the roles of a biogeographic barrier and environmental factors. J Biogeogr. 2022:49 (12 ):2296–2308. 10.1111/jbi.14504.
Milano  I, Babbucci  M, Cariani  A, Atanassova  M, Bekkevold  D, Carvalho  GR, Espiñeira  M, Fiorentino  F, Garofalo  G, Geffen  AJ, et al  Outlier SNP markers reveal fine-scale genetic structuring across European hake populations (Merluccius merluccius). Mol Ecol. 2014:23 (1 ):118–135. 10.1111/mec.12568.24138219
Miller  AD, Hoffmann  AA, Tan  MH, Young  M, Ahrens  C, Cocomazzo  M, Rattray  A, Ierodiaconou  DA, Treml  E, Sherman  CD. Local and regional scale habitat heterogeneity contribute to genetic adaptation in a commercially important marine mollusc (Haliotis rubra) from southeastern Australia. Mol Ecol. 2019:28 (12 ):3053–3072. 10.1111/mec.15128.31077479
Oigman-Pszczol  SS, Creed  JC. Quantification and classification of marine litter on beaches along Armação dos Búzios, Rio de Janeiro, Brazil. J Coast Res. 2007:232 :421–428. 10.2112/1551-5036(2007)23[421:QACOML]2.0.CO;2.
Oksanen  J, Blanchet  FG, Kindt  R, Legendre  P, Minchin  PR, O’hara  RB, Simpson  GL, Solymos  P, Henry  M, Stevens  H, et al  Package ‘vegan’.  Community Ecology Package. 2013:2 (9 ):1–295.
Pespeni  MH, Garfield  DA, Manier  MK, Palumbi  SR. Genome-wide polymorphisms show unexpected targets of natural selection. Proc Biol Sci.  2012:279 (1732 ):1412–1420. 10.1098/rspb.2011.1823.21993504
Pritchard  JK, Stephens  M, Donnelly  P. Inference of population structure using multilocus genotype data. Genetics. 2000:155 (2 ):945–959. 10.1093/genetics/155.2.945.10835412
Puillandre  N, Watkins  M, Olivera  BM. Evolution of Conus peptide genes: duplication and positive selection in the A-superfamily. J Mol Evol. 2010:70 (2 ):190–202. 10.1007/s00239-010-9321-7.20143226
Purcell  S, Neal  B, Todd-Brow  K, Thoma  L, Ferreir  MA, Bender  D, Maller  J, Sklar  P, de Bakker  PI, Daly  MJ, et al  PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007:81 (3 ):559–575. 10.1086/519795.17701901
Qiao  X, Hou  L, Wang  J, Jin  Y, Kong  N, Li  J, Wang  S, Wang  L, Song  L. Identification and characterization of an apoptosis-inducing factor 1 involved in apoptosis and immune defense of oyster, Crassostrea gigas. Fish Shellfish Immunol. 2021:119 :173–181. 10.1016/j.fsi.2021.09.016.34610453
Rajan  KC, Vengatesen  T. Molecular adaptation of molluscan biomineralisation to high-CO2 oceans—the known and the unknown. Mar Environ Res. 2020:155 :104883. 10.1016/j.marenvres.2020.104883.32072987
R Core Team . R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2013.
Reid  D . The genus Littoraria Griffith and Pidgeon, 1834 (Gastropoda: Littorinidae) in the Tropical Eastern Pacific. Veliger. 1999:42 :21–53.
Riesgo  A, Andrade  SC, Sharma  PP, Novo  M, Pérez-Porro  AR, Vahtera  V, González  VL, Kawauchi  GY, Giribet  G. Comparative description of ten transcriptomes of newly sequenced invertebrates and efficiency estimation of genomic sampling in non-model taxa. Front Zool. 2012:9 (1 ):33. 10.1186/1742-9994-9-33.23190771
Riginos  C, Crandall  ED, Liggins  L, Bongaerts  P, Treml  EA. Navigating the currents of seascape genomics: how spatial analyses can augment population genomic studies. Curr Zool. 2016:62 (6 ):581–601. 10.1093/cz/zow067.29491947
Romero  A, Novoa  B, Figueras  A. The complexity of apoptotic cell death in mollusks: an update. Fish Shellfish Immunol. 2015:46 (1 ):79–87. 10.1016/j.fsi.2015.03.038.25862972
Russo  J, Madec  L. Haemocyte apoptosis as a general cellular immune response of the snail, Lymnaea stagnalis, to a toxicant. Cell Tissue Res. 2007:328 (2 ):431–441. 10.1007/s00441-006-0353-7.17252246
San  L-Z, Liu  B-S, Liu  B, Zhu  K-C, Guo  L, Guo  H-Y, Zhang  N, Jiang  S-G, Zhang  D-C. Genome-wide association study reveals multiple novel SNPs and putative candidate genes associated with low oxygen tolerance in golden pompano Trachinotus ovatus (Linnaeus 1758). Aquaculture. 2021:544 :737098. 10.1016/j.aquaculture.2021.737098.
Santos  CA, Sonoda  GG, Cortez  T, Coutinho  LL, Andrade  SCS. Transcriptome expression of biomineralization genes in Littoraria flava gastropod in Brazilian rocky shore reveals evidence of local adaptation. Genome Biol Evol. 2021:13 (4 ):evab050. 10.1093/gbe/evab050.33720344
Sauerland  V, Kriest  I, Oschlies  A, Srivastav  A. Multiobjective calibration of a global biogeochemical ocean model against nutrients, oxygen, and oxygen minimum zones. J Adv Model Earth Sy. 2019:11 (5 ):1285–1308. 10.1029/2018MS001510.
Schaal  G, Riera  P, Leroux  C. Microscale variations of food web functioning within a rocky shore invertebrate community. Mar Biol. 2011:158 (3 ):623–630. 10.1007/s00227-010-1586-2.
Segovia  NI, González-Wevar  CA, Haye  PA. Signatures of local adaptation in the spatial genetic structure of the ascidian Pyura chilensis along the southeast Pacific coast. Sci Rep. 2020:10 (1 ):1–14. 10.1038/s41598-020-70798-1.31913322
Shapiro  MD, Marks  ME, Peichel  CL, Blackman  BK, Nereng  KS, Jónsson  B, Schluter  D, Kingsley  DM. Genetic and developmental basis of evolutionary pelvic reduction in threespine sticklebacks. Nature. 2004:428 (6984 ):717–723. 10.1038/nature02415.15085123
She  Z, Li  L, Meng  J, Jia  Z, Que  H, Zhang  G. Population resequencing reveals candidate genes associated with salinity adaptation of the Pacific oyster Crassostrea gigas. Sci Rep. 2018:8 (1 ):8683. 10.1038/s41598-018-26953-w.29875442
Shi  H, Wen  Z, Paull  D, Guo  M. A framework for quantifying the thermal buffering effect of microhabitats. Biol Conserv. 2016:204 :175–180. 10.1016/j.biocon.2016.11.006.
Sousa  R, Delgado  J, González  JA, Freitas  M, Henriques  P. Marine snails of the genus Phorcus: biology and ecology of sentinel species for human impacts on the rocky shores. Biol Res Water. 2018:2018 :141–147. 10.5772/intechopen.71614.
Sterza  JM, Fernandes  LL. Zooplankton community of the Vitória Bay Estuarine system (Southeastern Brazil): characterization during a three-year study. Braz J Oceanogr. 2006:54 (2-3 ):95–105. 10.1590/S1679-87592006000200001.
Storfer  A, Patton  A, Fraik  AK. Navigating the interface between landscape genetics and landscape genomics. Front Genet. 2018:9 :68. 10.3389/fgene.2018.00068.29593776
Takeuchi  T, Masaoka  T, Aoki  H, Koyanagi  R, Fujie  M, Satoh  N. Divergent northern and southern populations and demographic history of the pearl oyster in the western Pacific revealed with genomic SNPs. Evol Appl. 2020:13 (4 ):837–853. 10.1111/eva.12905.32211071
Tepolt  CK, Palumbi  SR. Transcriptome sequencing reveals both neutral and adaptive genome dynamics in a marine invader. Mol Ecol. 2015:24 (16 ):4145–4158. 10.1111/mec.13294.26118396
Tyberghein  L, Verbruggen  H, Pauly  K, Troupin  C, Mineur  F, De Clerck  O. Bio-ORACLE: a global environmental dataset for marine species distribution modelling. Global Ecol Biogeogr. 2012:21 (2 ):272–281. 10.1111/j.1466-8238.2011.00656.x.
Waldman  YY, Tuller  T, Keinan  A, Ruppin  E. Selection for translation efficiency on synonymous polymorphisms in recent human evolution. Genome Biol Evol. 2011:3 :749–761. 10.1093/gbe/evr076.21803767
Weir  BS, Cockerham  CC. Estimating F-statistics for the analysis of population structure. Evolution. 1984:38 (6 ):1358–1370. 10.1111/j.1558-5646.1984.tb05657.x.28563791
Westram  AM, Galindo  J, Alm Rosenblad  M, Grahame  JW, Panova  M, Butlin  RK. Do the same genes underlie parallel phenotypic divergence in different Littorina saxatilis populations?  Mol Ecol. 2014:23 (18 ):4603–4616. 10.1111/mec.12883.25113130
Yu  H, Li  Q. Genome-wide scan for positively selected genes in sessile molluscs highlights the genetic basis for their adaptation to attached lifestyle. J Ocean Univ China. 2018:17 (4 ):920–924. 10.1007/s11802-018-3611-x.
Zhbannikov  IY, Hunter  SS, Foster  JA, Settles  ML.  2017. SeqyClean: a pipeline for high-throughput sequence data preprocessing. In Proceedings of the 8th ACM International Conference on Bioinformatics, Computational Biology, and Health Informatics: 407-416.
Zhong  X, Li  Q, Yu  H, Kong  L. SNP mining in Crassostrea gigas EST data: transferability to four other Crassostrea species, phylogenetic inferences and outlier SNPs under selection. PLoS One. 2014:9 (9 ):e108256. 10.1371/journal.pone.0108256.25238392
