Outside your shell: exploring genetic variation in two congeneric marine snails across a latitude and temperature gradient

Abstract

Intertidal organisms withstand extreme temperature fluctuations, and theirability to cope with this variation may affect their distributions across the seascape. Genetic variation and local environments likely interact to determine variation in thermal performances across intertidal species’ ranges, so characterizing the relationship between temperature variation and population structure is key to understanding the biology of marine invertebrates. Here, we use 2bRAD-sequencing to examine population genetic structure in two congeneric intertidal marine gastropods (Crepidula fornicata, C. plana), sampled from locations along a natural temperature gradient on the Northeast shores of the United States. These two species share similar life histories, yet C. plana exhibits a narrower distribution than C. fornicata. Our results demonstrate that both species show patterns of genetic divergence consistent with isolation by distance, though this pattern was only significant in C. fornicata. Both putatively selected and neutral loci displayed significant spatial structuring in C. fornicata; however, only putatively selected loci showed significant clustering in C. plana. When exploring whether temperature differences explained genetic differentiation, we found that 9–12% of genetic differentiation was explained by temperature variation in each species even when controlling for latitude and neutral population structure. Our results suggest that temperature shapes adaptive variation across the seascape in both Crepidula species and encourages further research to differentiate our results from models of neutral evolutionary drift.

Introduction

Species distributions are governed by a variety of biotic and abiotic factors. In marine taxa, many species disperse through a bipartite lifecycle, which includes a sessile adult stage and a planktonic larval stage. Historically, it had been assumed that these highly dispersive larvae lead to panmixia between populations, even across large spatial scales [1, 2]. But advancements in modern genomic sequencing technologies for non-model systems have challenged this view and research has demonstrated that many marine organisms are much more genetically structured than previously appreciated (for review, see: [3]). This genetic structure can manifest through neutral or selective processes. Neutral processes can structure genetic variation through differences in migration rates [4] as well as through oceanographic properties that facilitate or constrain larval dispersal across the seascape [5]. However, even in the face of high gene flow, selection can also shape allele frequencies following environmental clines [6]. Disentangling how these evolutionary processes operate across environmental gradients is key baseline knowledge for understanding how environments shape species distributions and how these distributions might shift under climate change.

Temperature is an important abiotic factor shaping species distributions as temperatures constrain many physiological processes [7]. Due to its ubiquitous effects, temperature leaves broad fingerprints on the genetic architecture of species and populations [8]. Temperature can act as a selective agent due to its ability to govern the basic functions of numerous cellular processes [9], which can ultimately shape patterns of local adaptation [10]. Perhaps less intuitively, temperature regimes can theoretically drive population genetic patterns via neutral processes. Large declines in effective population sizes can lead to drift becoming a stronger evolutionary force relative to selection [11], which could be a product of a major heatwave [12]. As climate change continues to warm habitats globally [13], these warmer temperatures will be a critical factor determining genomic variation and structure both within and across species in the coming decades.

Temperatures naturally vary along latitudinal gradients and this variation can be leveraged to investigate how populations will respond to temperature changes [14]. For example, latitudinal gradients can be used to generate and inform models aiming to predict species range expansion limits [15], or they can be used to isolate the genomic determinants of thermal tolerance [16, 17]. Signs of local adaptation can follow predictable environmental gradients, even in marine organisms with high dispersal potential (for review see [18]). For example, in a reciprocal crossing experiment in the eastern oyster (Crassostrea virginica), evidence of local adaptation between northern and southern genetic lineages across the eastern Florida coast was found [19]. This pattern suggests that natural selection maintains genetic step-cline patterns following a latitude gradient despite high dispersal potential in this species. In addition, due to the highly predictable nature of temperature differences along latitudinal gradients, these clines may be key for identifying ecologically-relevant loci for thermal adaptation [20]. However, it is important to note that thermally adaptive loci may lead to fixed phenotypes or plastic phenotypes with underlying thermally adapted gene expression [21]. Phenotypic plasticity also enables organisms of the same genotype to occupy divergent thermal regimes [22], and this trait is thought to be particularly relevant in intertidal species with wide spatial ranges that experience high environmental variation within their habitats and across their range. These latitudinal gradients may therefore provide important insights into the ecological and evolutionary processes shaping intertidal species distributions.

Many studies of population structure in marine invertebrates have historically used a small number of markers, such as microsatellites or mitochondrial DNA [23]. While these loci can inform coarse scale relationships between temperature and genetic structure [24, 25], using a limited number of markers may provide an incomplete picture of more fine scale patterns. Previous work has often detected marine invertebrate population structure across wide spatial scales of hundreds to thousands of kilometers (for review see [26]), but the low resolution provided by only a few loci makes it challenging to relate this variation to environmental parameters. For example, Schmidt et al. [27] was unable to detect genetic subdivision among sampled populations of the marine snail Littorina obtusata using microsatellite markers despite strong evidence for local adaptation along a latitudinal gradient. However, when they included a known a priori candidate allozyme marker they detected significant differentiation of this marker following the latitudinal cline [27]. This challenge may be amplified in species with high dispersal capabilities and potentially high gene flow (as in the case for many marine invertebrates), as this can restrict the signal of adaptive loci to just a few important regions. Genomic studies leveraging more modern technologies can provide finer resolution [28], allowing for more direct assessments of the relationship between environmental and genomic variation.

Much like how latitude can be leveraged to compare thermal regimes, comparing congeners offers the opportunity to understand conserved mechanisms involved in response to environments. For example, a comparison of congeneric porcelain crabs across a thermal gradient established the role of different morphological and biochemical features that lead to thermal tolerance [29]. Furthermore, contrasting congeners can elucidate links between thermal tolerance and species distributions. For example, by comparing two invasive and two non-invasive beach grasses, it was found that the two species with broad distributions and high invasion success exhibited wider thermal breadths than the other two species with more narrow ranges and low invasion success [30]. Finally, population genetic comparisons of congeners offer important insights into how genetic variation and connectivity are shaped by seascapes [31] and because closely related species often share similar substitution and mutation rates [32], this allows for more direct comparisons of gene flow.

Two congeneric species of Crepidula (C. fornicata, C. plana) live in sympatry along the northeast Atlantic coast of the United States and offer a promising comparative system to characterize how temperature shapes genetic variation through neutral and selective forces. These sister species of benthic marine gastropods are commonly known as ‘slipper snails’ and have a protandric hermaphrodite lifestyle with smaller males forming stacks of varying densities on top of larger females [33]. While both species exhibit stacking behaviour, only C. fornicata forms large gregarious stacks. Furthermore, C. plana are smaller and often found on the undersides of rocks or shells and are commonly found in symbiosis with species of hermit crabs [34]. Recent attention has focused on C. fornicata as an emerging model system for Lophotrochozoan development [35], as well its use in CRISPR/Cas-9 mediated gene knockouts [36]. While a promising model system, C. fornicata has received most of its attention due to its invasive success in Europe and on the west coast of North America [37], which represents a key difference between its native congener. Despite sharing similar life history characteristics such as brooding larval development [38, 39] and recorded opportunities for invasion in Europe [40] C. plana has not successfully invaded Europe and has a more restricted northern range limit relative to C. fornicata [41]. Invasive success in C. plana is limited, despite sharing similar life history characteristics and recorded opportunities for invasion in Europe [40]. While C. plana was originally thought to have the same spatial distribution as C. fornicata [42], sequencing has uncovered three distinct cryptic species within the C. plana complex: C. plana, C. depressa and C. atrasolea [43]. These cryptic species of C. plana separate between known biogeographic breaks where current dynamics form natural dispersal barriers for many marine organisms [5] and hereafter any reference to C. plana refers to the species within the complex, not the complex itself. It is possible that C. plana’s pelagic larval duration (PLD) contributes to its truncated range given that its maximum PLD is lower than C. fornicata’s (12 days vs. 21 days; [44]. The duration of an organism’s PLD can also contribute to genetic structure, where shorter PLD leads to greater differentiation [45]. However, a meta-analysis on reef invertebrate connectivity found little support for PLD correlating with genetic structure [23]. It remains unclear whether PLD differences between congeneric Crepidula explain variation in range sizes or whether other biotic or abiotic factors explain these differences.

Here we compare how genetic variation between the congeneric gastropods C. fornicata and C. plana is partitioned across five sites along the northern portion of their native range in the United States using reduced representation DNA sequencing (2bRAD-seq). We then explored how temperature variation may explain patterns of genetic diversity across this latitudinal gradient. Our findings contribute to growing research exploring how different evolutionary processes drive genetic variation in marine organisms across latitudinal gradients.

Methods

Sample collection and environmental data

During the low tides of July and August 2021, individuals identified as Crepidula fornicata and C. plana were collected from five different sites along the Northeast shores of the United States (Fig. 1A). Animals were removed from intertidal substrate using a flathead screwdriver and transported in a cooler containing continuously aerated sea water. No C. plana were found at the northernmost sampling site (Robbinston, ME). All animals were transferred into holding aquaria at Boston University’s Marine Invertebrate Research Facility within 24 h of collection, before being removed from their shells with tweezers, preserved in 200 proof ethanol and stored at -80˚C until processing. To characterize seawater and air temperature profiles across sites, hourly temperature data between January 2020 to January 2021 were obtained from the National Oceanic and Atmospheric Administration (NOAA) weather buoys that were closest to each sample location (Table S1). These temperature data were plotted along with a local regression line (loess smoothing) by site location to illustrate seasonal variation in temperature across sites (Fig. 1B-C). Average temperatures for each site were calculated for spring (March, April, May), summer (June, July, August), fall (September, October, November), and winter (December, January, February) and these values were used for downstream association analyses.

Fig. 1

Collection sites and associated temperature profiles for Crepidula fornicata and C. plana. A) Sampling locations of Crepidula spp. along the northeast shore of the United States. Mean hourly temperature plotted along with bolded local regression line (loess smoothing) collected from NOAA weather buoys (Table A3-1) from January 2020 - January 2021 for B) sea surface temperature, and C) air temperature

2bRADseq library preparation

A small portion of the animal’s foot (approximately the size of a rice grain) was removed with a scalpel and DNA was isolated using Omega BioTek EZNA Mollusc DNA kit following manufacturer’s instructions with one additional DNA washing step (three washing steps total). DNA was normalized to 100 ng in 4uL and prepared for 2b-RAD sequencing following [46]. Six samples were prepared in duplicate to serve as technical replicates to assist with downstream bioinformatic analyses. A total of 173 animals (see Table S1 for sample sizes across species and sites) were barcoded and sequenced across three lanes of Illumina HiSeq 2500 using single end 50 bp at Tufts University Core Facility.

Filtering, quality control and single nucleotide polymorphism identification

Per-sample library sizes and de novo mapping rates are available in supplemental Table S2. Raw reads were deduplicated, trimmed and quality filtered using FASTX toolkit (http://hannonlab.cshl.edu/fastx_toolkit) with Phred scores > 15 (-q 15,15 -m 36). Since no reference genomes are publicly available for either C. plana or C. fornicata, de novo references were created for each species separately. To construct these references, common tags (major alleles) were collected with the minInd filter set as 10% of sample numbers and then any tag having more than seven observations without a reverse-complement was discarded. References were made using bowtie2 [47] setting the chromosome number to 17, based on known chromosome numbers for a related species Crepidula unguiformis [48]. Sequence alignment files (SAM) from all samples, including technical replicates, were then compressed into BAM files using samtools (version 1.12). Given that depth of coverage for both species was < 10X, we utilized angsd (version 0.935; for filters used see Table S3) which is designed to analyze low depth of coverage sequence data to identify single nucleotide polymorphisms (SNPs) by inferring genotype likelihoods into downstream analyses [49]. Seven of 173 samples (4% of samples; Table S4) had > 85% missing data and were subsequently discarded. Angsd was re-run after these low-coverage samples were removed (for filters see Table S4). The clustering of technical replicates was verified on a dendrogram with hierarchical clustering using the identity by state matrix output from angsd using hclust in R (version 4.0.2; Figure S1). All technical replicates showed strong clustering (< 0.2 height) in both C. fornicata and C. plana independent references with initial angsd minInd filter set to 80% of respective sample sizes. Additionally, a subsample of variant sites was taken from C. plana to match the total SNP number for C. fornicata. Neutral, and outlier SNPs were identified, and population structure was determined from this downsampled list of loci using the same pipeline and filters as described below.

Identifying loci under selection

Candidate loci that are putatively differentiated within a species due to selection were determined using pcadapt [50]. This method identifies candidate SNPs as outliers with respect to population structure. This structure was therefore first inferred using a principal component analysis (PCA) using PCAngsd (version 1.11; [51]), which estimated a covariance matrix of individual allele frequencies from the genotype likelihood files generated by angsd described above. PCAngsd also generated principal component selection statistics, which were used by pcadapt to perform a genome-wide selection scan to identify outlier SNPs. Outliers for each dataset were determined by the Mahalanobis distance of z-scores along the first principal component followed by a χ² test with the outlier SNP alpha set as p < 0.001 (i.e., significantly different from mean distance for each locus in PCA space). Neutral SNPs were identified using this same method, with a SNP alpha set to p > 0.05. While more conservative thresholds (e.g., p > 0.20) would increase confidence in neutrality of these SNPs, we chose this cutoff to retain a broader set of putatively neutral loci for downstream analysis. The lists of outlier and neutral SNPs for all datasets were then separately re-run in angsd using the same filters as above to determine how genetic variance is shaped when exploring (i) all loci, (ii) neutral loci, or (iii) loci under selection.

Assessing population structure, genetic diversity and demographic history

Within each species a multidimensional scaling (MDS) analysis was used plotting the Euclidean distances between all samples. This analysis was performed on a covariance matrix for (i) all loci, (ii) neutral loci and (iii) outlier loci based on single-read resampling calculated in angsd described above. The position of samples along the first and second MDS axes were visualized, and a permutation multivariate analysis of variance determined significant differences in clustering based on sample site and species using the adonis function in vegan [52]. Admixture between individuals within a species was explored using NGSadmix [53], which analyzes genotype likelihood files from all SNPs and neutral SNPs only. Optimal K values were determined using the Evanno method [54] by bootstrapping the above analyses 10 times for each K value between 2 and 5, and using the log likelihood of the estimates from each iteration as input into CLUMPAK [55]. To determine genetic differentiation within each species between collection sites, angsd calculated site allele frequencies (SAFs) for each population which were visualized by plotting histograms of the SAFs between species and shared sampling locations. Then, realSFS calculated the folded 2-dimensional site frequency spectrums (SFSs) for all possible pairwise comparisons. These folded SFSs were used to calculate Global Watterson’s θ, π, and Tajima’s D for each population using the thetaStat tool in angsd. Mean differences in these metrics among sampling locations and species were assessed using a Dunn’s test with multiple test correction using the Bonferroni method [56]. SFSs were then used as priors with the SAF to calculate pairwise F_ST. Only weighted pairwise F_ST values between populations within a species are reported here. To test for patterns of isolation by distance (IBD), a matrix of geographic distances (shortest straight line) between sites was correlated with a matrix of pairwise F_ST values using a Mantel’s test with all possible permutations implemented in vegan. Finally, to reconstruct the history of effective population size changes on our folded SFS, Stairway Plot v2 [57] was used, which is an unsupervised method requiring no prespecified demographic models. For all demographic analyses we used default parameters for mutation rate (1.2e-8 per base per year), and a generation time of 3 years.

Identifying loci associated with temperature regimes

Ordinances from the MDS analyses for each species were correlated with the average seasonal temperature data obtained from NOAA weather buoys (see Table S1). Locations on the first axis of each MDS, which explained the most genetic variance, were plotted against temperature regimes using linear regressions. A series of redundancy analyses (RDAs) were then used to model the relationships between genomic variation with environmental predictors and latitude. These were performed using genotype matrices for all loci and neutral loci and in cases of missing data the most common genotype at each SNP was assigned. The genotype matrix of all loci was used as the response variable and the neutral genotype matrix was used to assess the contribution of neutral population structure in the RDAs. To assign neutral population structure, we first performed a principal component analysis on neutral loci and then used eigenvalues from the first principal component as a variable in our models. Significant environmental variables were chosen based on a forward selection strategy with the ordi2step function in vegan [52] with a significance value of p < 0.01 and 1000 permutations. A full model was chosen where genetic variation was predicted from significant environmental variables whose model from the above step had the lowest Akaike information criteria (AIC). This model then further included latitude and neutral genetic structure and was compared with a pure climate model (significant environmental variables as explanatory and latitude and neutral population structure as conditional), pure latitude model (latitude as explanatory and environmental and neutral population structure as conditional), and a pure neutral population structure model (neutral population structure as explanatory with environmental and latitude as conditional) using partial RDAs.

Results

Across all samples, sequencing resulted in an average of 693,161 (SE = 28,787) reads per individual, which was reduced to 617,260 (SE = 26,182) after filtering and mapping to species specific de novo references. Within species, mapping of C. fornicata (reference: 14.16 Mbp; 7.08X coverage) identified 4,379 total SNPs (outliers: 225; neutral: 3,570). C. plana (reference: 15.81 Mbp; 6.85X coverage) had approximately 4.5x as many SNPs (outliers: 317; neutral:17,322).

Population structure

Multidimensional scaling (MDS) plots demonstrate significant population structure across sites for C. fornicata regardless of the loci being tested (All = 4,379, Neutral = 3,570, Outlier = 225) (Fig. 2; A-C, p_site < 0.001). In contrast, C. plana exhibited significant population structure for all loci (p_site < 0.001) and outlier loci (p_site < 0.001), but no structure was observed for neutral loci (p_site = 0.771) (Fig. 2; D-F). Population structure was further explored using all loci and neutral loci using ADMIXTURE, which identified (K = 3) ancestral populations in both C. fornicata (Figure S2) and C. plana (Figure S3). These admixture analyses showcased ancestral population assignments shifting along the latitudinal gradient (Figure S2, S3). When pairwise genetic divergences were calculated between sites (Table S5; F_ST), all comparisons were significant, but low in magnitude (< 0.018), with the largest F_ST being between sites with the greatest distance in both species. This pattern resulted in significant isolation by distance (IBD) for C. fornicata (Fig. 3; Mantel’s r = 0.760, p = 0.0167), but not for C. plana (Mantel’s r = 0.658, p = 0.083). However, this IBD signal in C. fornicata was no longer significant when the most northern sampling site (Robbinston, ME) was removed from the C. fornicata analysis (r = 0.627, p = 0.083) and both the slope and correlation-coefficient were similar in both species, suggesting that the C. plana dataset was potentially underpowered. Furthermore, downsampling the full panel of SNPs to 4,379 random variant sites in C. plana led to 87% of SNPs being neutral (3082 SNPs) and 1.9% identified as outliers (85 SNPs). These were similar proportions to the full data set (19,820 variant sites), where 87% were identified as neutral and 1.6% as outliers. However, our ability to detect genetic structure differed between the downsampled and all variant datasets as we found significant structuring in neutral but not outlier loci in the downsampled dataset (Figure S4), which was the opposite pattern when the full SNP dataset was used (Fig. 2).

Fig. 2

Multidimensional scaling (MDS) plots showing population genetic structure of Crepidula fornicata (A-C) and C. plana (D-F). Clustering is shown for single nucleotide polymorphisms (SNPs) across all loci (A, D), and those that were identified as neutral (p > 0.05; B, E), and outlier (p < 0.001; C, F) SNPs via pcadapt. Significant genetic structure was detected across sampling sites using a permutational multivariate analysis of variance in all analyses except for neutral loci in C. plana (denoted with grey box)

Fig. 3

Estimated isolation by distance (IBD) as measured by the relationship between pairwise physical distances between sites (km) and pairwise genetic divergences (F_ST) in Crepidula fornicata (purple) and C. plana (pink). IBD was calculated using a Mantel’s test. C. fornicata exhibited significant IBD, but C. plana did not, potentially due to its relatively lower power (see text)

Genetic diversity

Global Watterson’s θ was greater in C. plana (Kruskal-Wallis χ² = 5.54, df = 1, p = 0.0186); however, no difference in genetic diversity (π) between species was observed (Kruskal-Wallis χ² = 2.42, df = 1, p = 0.120). C. fornicata had a greater average Tajima’s D across all populations (Kruskal-Wallis χ² = 14.0, df = 1, p < 0.001). Comparisons of the sampling locations within species show no differences in π between locations in both C. fornicata and C. plana. Within C. fornicata, Watterson’s θ was significantly lower and Tajima’s D higher in Robbinston, ME compared with Newport, RI and Cape May, NJ (Fig. 4B-C). Within C. plana, Watterson’s θ was significantly lower and Tajima’s D higher in Beverly, ME compared to all other sampling locations (Fig. 4B-C).

Fig. 4

Genetic diversity estimates. A) genetic diversity (π), B) Watterson’s θ, and C) Tajima’s D. Points represent means, with error bars denoting standard error. Significant pairwise comparisons (p_adj < 0.05) within species are denoted by dashed horizontal lines for Crepidula fornicata and solid lines for C. plana

Demographic inference

Stairway plot analysis of effective population (N_e) sizes of C. fornicata suggest recent N_e contractions over the past few thousand generations at all sites except for Newport, RI (Fig. 5A). These contractions vary in their magnitude and follow the latitudinal gradient with Robbinston, ME having the smallest contemporary N_e and increasing at each site moving southward. This contrasts with C. plana’s results, which suggest increasing N_e over 100,000 generations ago followed by relative stasis, except for Beverly, MA (Fig. 5B). In this location, the trend follows a similar pattern to C. fornicata with contemporary contractions in N_e.

Fig. 5

Stairway plot of inferred effective population (N_e) sizes over time for A) Crepidula fornicata and B) C. plana. Lines show median N_e and ribbons represent 95% confidence intervals. Note that y-axis scales differ between species

Associations of temperature and genetic variance

Two approaches were used to measure the association of temperature with genetic variation. The first approach associated ordinations along the first MDS axis for (i) neutral, and (ii) outlier loci with average seasonal (spring, summer, winter, fall) air and water temperatures. Using this method, structure in neutral loci was correlated with seasonal temperatures in C. fornicata (p < 0.05), except for winter water temperatures not showing a significant association (Figure S5). This pattern contrasted with C. plana where structure in neutral loci exhibited no associations with seasonal water and air temperatures, except for summer and spring air temperatures exhibiting weak associations (Figure S6). Lastly, structure in outlier loci showcased the strongest associations with water and air temperatures in both species. In C. fornicata, structure observed in outlier loci was correlated with temperatures across seasons (Figure S7; all p < 0.001) and similar patterns were observed in C. plana except for spring sea surface temperature having no significant association (Figure S8).

While the above approach assesses correlations between genetic variation and temperature, we also employed a multivariate approach to further characterize how temperature and spatial heterogeneity impact genetic variation using a series of redundancy analyses (RDAs; Fig. 6). When using a forward selection strategy to select a pure temperature model that best explained our data, we found that summer sea surface temperature in C. fornicata and fall sea surface temperatures in C. plana were the best fit models (lowest AIC). When these top environmental models were used as the explanatory variables, while controlling for latitude and neutral population structure in an RDA, temperature contributed 9% of the explainable variance in C. fornicata and 12% in C. plana. Latitude was also a significant predictor of the explainable variance in both species (10% C. fornicata, 12% C. plana) when controlling for temperature and neutral population structure. Neutral population structure explained the most variation in both species (C. fornicata: p < 0.05, 77%; C. plana: p < 0.05, 71%), when controlling for latitude.

Fig. 6

Environmental variation associated with genetic variation in all loci of A) Crepidula fornicata and B) C. plana. Bar plot denotes variance explained in redundancy analysis (RDA). Confounded percentage denotes variation that could not be separated between temperature models (summer temperature in C. fornicata and fall temperature in C. plana) and geographic models (latitude). C) Average summer water temperatures were a significant driver of variation for C. fornicata and were correlated with genetic variation of outlier loci using the first MDS axis from Fig. 3. D) Average fall water temperatures were a significant driver for C. plana genetic variation within outlier loci using the first MDS axis from Fig. 3

Discussion

We investigated patterns of genetic variation in the two congeneric marine gastropods Crepidula fornicata and C. plana across a latitudinal gradient (approximately 900 km) spanning the northern portion of their native ranges. These species share many similar life history traits, such as planktonic larval development, yet they exhibit key differences in the extent of their spatial ranges. Relative to C. plana, C. fornicata is found at higher and lower latitudes across its local range [43] and C. fornicata also exhibits contemporary invasive ranges globally [37]. We inferred demographic histories across this latitudinal gradient and tested for associations between genetic variation and temperature. Broadly, similar divergence patterns were observed between C. fornicata and C. plana across latitude, with C. fornicata exhibiting higher genetic divergence between locations. Our demographic analyses revealed that contemporary contractions of effective population (N_e) sizes in C. fornicata followed a latitudinal gradient where N_e contractions were higher at higher latitude sites. In contrast, C. plana population sizes have remained much more constant. However, in this analysis we assume the same mutation rate and given our polymorphism differences and potential for different mutation rates we highlight that this may present a bias in our interpretations. Finally, we find evidence that temperature regimes are associated with genetic variation in both neutral and outlier loci in C. fornicata, but only outlier loci in C. plana. These results suggest subtle differences in genetic differentiation between these congeners, perhaps due to differences in gene flow between populations or historic bottlenecks in C. fornicata that were not experienced by C. plana. Our results support a possible role for temperature in shaping divergence patterns across the seascape, though further work disentangling temperature from other abiotic factors that are correlated with latitude is needed.

Interestingly, C. plana had approximately 4.5X as many total SNPs as C. fornicata despite similar numbers of outlier SNPs (C. fornicata: 225, C. plana: 317), similar depth of coverage (C. fornicata: 7.08X, C. plana: 6.85X) and similar total numbers of invariant and variant loci (C. fornicata: 14.16Mbp, C. plana: 15.81Mbp). As the number of variant sites may impact our power to resolve structure, we explored a reduced dataset of C. plana by downsampling loci to be consistent with C. fornicata. Downsampling had an impact on our ability to detect genetic structure in both neutral and outlier sites (Figure S4). We were not able to determine the reason why these polymorphism differences between species exist, but we caution that results can be influenced by these power imbalances in comparative genomics. Nevertheless, it is reasonable to suspect biological explanations of these differences, which may include different mutation rates [58] or variation in mechanisms of rare allele maintenance or removal across the seascape [59]. We did observe an excess of rare alleles in C. plana relative to C. fornicata, which could contribute to more total SNPs detected in C. plana (Figures S10). Further research would be needed to confirm these polymorphism differences between species and to explore how these differences may impact the power to detect population structure.

We detected significant genome-wide population structure across the latitudinal gradient in all, putatively neutral, and selected loci in C. fornicata (Fig. 2). Similar patterns emerged in C. plana in all and selected loci; however, there was no significant structuring in neutral loci (Fig. 2). Downsampling the number of variant sites in C. plana led to a reversal in this pattern with population structure detected in neutral loci but not outlier loci (Figure S4). Further research is clearly needed to determine why there are differences in polymorphisms between these species and how this may impact population structure. Nevertheless, the patterns using our full data were consistent with our admixture results where a distinct north to south gradient in population structure for all loci and neutral loci was observed in C. fornicata (Figure S2) whereas C. plana only exhibited latitudinal associations in all loci and not neutral loci (Figure S3). While population structure was observed, low between-site F_ST was also found in both species, suggesting high gene flow across this range (Table S5). We found a significant pattern of isolation by distance (IBD; [60]), in C. fornicata but not in C. plana (Fig. 3); however both species exhibit similar IBD characteristics, and this trend could potentially become significant with additional sites sampled. While different structure patterns between species could be driven by a discrepancy of the total number of loci, we do not think this is the case here as we found similar numbers of total loci in the de novo references in both species. Overall, each of our analyses paints a consistent picture that population structure along the latitudinal gradient exists but is relatively weaker in C. plana.

These findings are largely consistent with previous work in Crepidula, which have explored population structure using few loci. For example, genetic structure in C. fornicata but not C. plana was found along the coast of New England using cytochrome oxidase I [43] and allozymes [61]. Our C. fornicata results are also similar to those observed using microsatellite markers across their native range [62, 63]. Together, our results provide further evidence that biophysical features across C. fornicata’s native range lead to genetic divergence across the seascape, and these same features impact C. plana, although to a lesser degree. Work in other intertidal invertebrates such as in mussels (Mytilus edulis; [64]), oysters (Crassostrea virginicacite; [65]), and barnacles (Pollicipes pollicipes; [66]) have all consistently shown evidence for neutral population structure. It is unclear if there are specific biological features of C. plana that facilitate higher connectivity and reduce neutral population structure. Alternatively, as we show with our downsampled analysis, the degree of polymorphism within a species could also impact these findings.

In addition to population structure, we also measured differences in genetic diversity as well as inferred demographic histories across the latitudinal gradient in both species. In C. fornicata we found lower Watterson’s θ and higher Tajima’s D in samples from the most northern sampling site relative to two southern locations (Newport, RI, Cape May, NJ; Fig. 4B-C). This pattern suggests a depletion of rare alleles in the north and may correspond with inferred decreases in N_e over the last few thousand years (Fig. 5A). This recent population contraction is contrary to demographic signals inferred in C. plana as well as other temperate marine invertebrates such as scallops whose northern populations show steady N_e expansions over the last few thousand years [67]. C. plana maintained consistently high N_e at all sites except Beverly, MA, which had a much smaller contemporary N_e and demographic history that more closely mirrors pattern observed in C. fornicata. This site is a clear outlier in both diversity statistics as well as in our demographic analysis. Our demographic analyses are consistent with a model in which C. plana populations have remained resident to the northeast, while C. fornicata populations either declined during the last glacial maximum or reinvaded as glaciers receded. These demographic differences are interesting given how these species share a similar biogeographical history and may hint to differences in their ecological or evolutionary history. For example, recent northern population contractions may lead to increased signals of IBD in C. fornicata due to founder effects not experienced by C. plana. In this model, the northern populations represent the leading edge of an expanding wavefront, perhaps beginning at the end of the last glacial maximum, rather than a contracting native population. Future research using annotated whole genome sequences and haplotype patterns is warranted to further clarify the specific drivers of differentiation in these populations over their evolutionary histories.

There are many plausible explanations for the differences observed here between C. fornicata and C. plana. C. plana was not observed at the northernmost site, perhaps because it has a smaller thermal niche width than C. fornicata. A more truncated niche space could lead to patterns of population structure in C. plana to become less pronounced simply because it has a smaller geographic range. Furthermore, C. plana is smaller and is often found in association with hermit crabs [34], which could increase their mobility and connectivity potential. In contrast, C. fornicata are larger and tend to exhibit limited movement with the exception of small males [39]. However, a key difference in reproductive behaviours between these species is that C. plana are more often found in pairs and tend not to form large, gregarious stacks like C. fornicata. These stacks function as independent mating groups, often with skewed sex ratios with more males than females [68]. Therefore, despite a larger census population size, C. fornicata’s stacking behaviour could be a mechanism to lower their Ne. This reproductive behavior could increase inbreeding potential and the influence of genetic drift, which would both increase neutral genetic divergence between populations [69].While this is certainly the case in gonochoric organisms [70], hermaphroditic organisms may be able to overcome these disadvantages as individuals can breed as both sexes. For example, work in fishes has demonstrated that sequential hermaphroditism is an evolutionary stable strategy that minimizes the evolutionary constraints (such as reduced N_e) of a skewed sex ratio [71]. Crepidula therefore presents opportunities for future research to understand how these unique reproductive characteristics may impact effective population size, rates of inbreeding, or other aspects of biology that may drive the population genetic patterns observed between C. fornicata and C. plana.

Genetic differentiation at putative loci under selection associates with temperature in C. fornicata and C. plana

Species in intertidal zones experience frequent fluctuations in temperature and environmental conditions due to the changing tides [72]. These fluctuations become more pronounced at higher latitudes, with both temperature and tidal extremes intensifying. Therefore, adaptations to broad thermal regimes are necessary to survive these harsh conditions [73]. Here, we find support for temperature shaping genetic variation across a latitudinal gradient in both Crepidula species, and this pattern was most prominent in putatively selected loci. Correlations between temperature and MDS1 of neutral loci were either non-significant or weakly significant in both C. fornicata (Figure S5) and C. plana (Figure S6). In contrast, associations using outlier (putatively selected; Figure S7, S8) loci alone were considerably stronger in both species. Therefore, these loci might be driven by adaptation to thermal regimes rather than through neutral processes, though we cannot rule out neutral processes completely (see discussion of RDA below). Previous work by Thieltges et al., (2004) concluded that winter temperatures were the primary constraint for range expansion in northern populations of C. fornicata [74]. From their results, we might expect that winter months would impart strong selective pressures and therefore be strongly correlated with outlier loci observed here. It has also been proposed that cold water acts as a significant selective force in C. plana, with limitations in cold water tolerance potentially explaining why C. plana has a more limited northern range and has failed as an invasive species in Europe [75]. Our results in both analyses fail to support this prediction as the association between temperature and differentiation amongst outlier loci was strongest during the summer (C. fornicata; Fig. 6A) and fall (C. plana; Fig. 6B) when water temperatures are the warmest in this region. Together, our data suggest that patterns of adaptive genetic divergence in Crepidula are more strongly associated with warmer temperatures. This could indicate thermal adaptations or be linked to reproduction or larval development as spawning occurs in the warmer months in both species [44, 76]. Future research characterizing the thermal performance curves of these populations and across life-stages coupled with high resolution whole genome sequencing could provide further resolution into the biological mechanisms impacted by temperature in these species.

Genomic signatures of selection can sometimes be falsely attributed to variance that is evolutionarily neutral [77]. A limitation of the correlation analyses we performed between outlier loci and temperature variables is that outlier loci are detected because they have larger than expected frequency differences between populations, which means their latitudinal signals are inherently stronger. Therefore, we controlled for neutral population structure by performing a redundancy analysis, which is a more conservative approach to identify signals of selection [78]. We found that after controlling for latitude and neutral population structure, temperature explained approximately 9% of explainable variance in C. fornicata and 12% in C. plana (Fig. 6). While temperature explained similar amounts of genetic variance in both species, our redundancy analysis highlights that different seasons contribute to this variance between congeners. Summer (C. fornicata) and fall (C. plana) water temperatures explained the most variation, suggesting that warmer waters are stronger drivers of this genetic variation. This pattern contrasts with many terrestrial species where winter lows have the largest explanatory power [79]. While this may be true for terrestrial plants, likely due to frost and freezing damage, this aspect of temperature may not be as relevant in marine organisms. Indeed, RDA approaches often identify warm water as explaining adaptive genetic variance in marine invertebrates at latitudes like those sampled here. For example, warmer sea temperatures explained genetic variation in sea cucumbers (Parastichopus californicus; [80]), and urchins (Strongylocentrotus purpuratus; [81]) along latitudinal clines. However, we only found that 9–12% of the explainable variance was associated with temperature. A more exhaustive model that includes accurate daily temperature and tidal variations would better inform the RDA analysis. Here, we used oceanic weather buoy data, which is a coarse estimate of the temperature variation experienced by these intertidal organisms. Future work should consider capturing more fine scale data from the intertidal environment of each location. For example, deployment of biomimetic “robosnails” has been used amongst C. fornicata populations in Rhode Island, USA and better capture extreme diel temperature variability (20–43˚C; [73]). Additionally the use of satellite-based weather data may better characterize differences in the abiotic environment between sampling locations and could be used to better inform RDA models. Nevertheless, despite these limitations we show evidence for temperature explaining genetic variance and this variation is likely adaptive.

Previous work along the northeast coast of the United States and Canada have found evidence of structured genetic variation along this latitudinal cline. For example, sea scallops (Placopecten magellanicus; [82]) and the marine snail (Macoma petalum; [83]) both demonstrate distinct north and south population structure separated by the Bay of Fundy, Canada and this structure was most pronounced in loci putatively under selection [84]. Further exploration of sea scallop populations across warmer inshore and cooler offshore populations, found that temperature differences between environments explained genetic structure [85]. Overall, our data support the consensus that population structure of marine invertebrates in this region may be influenced by varying temperatures.

Conclusion

We compared population genetic patterns in two congeneric marine gastropods with similar life history characteristics in the genus Crepidula to explore how genetic variation is partitioned across a wide latitudinal gradient. We found parallel patterns of genetic divergence and diversity between species, with subtle differences in connectivity across the seascape. We therefore predict that under future ocean warming associated with climate change, these warmer temperatures will continue to influence genetic variation in these species.