DISCUSSION
We implemented a seascape genomics approach to test for differences in neutral and adaptive genetic connectivity between two intertidal marine gastropods, the direct developer Nucella lapillus , and the broadcast spawning Steromphala umbilicali s, within similarly distributed ranges in the UK and Ireland. Neutral genetic structure conformed to expectations of lower connectivity in the direct developing species, with pairwise FST being on average 11.4x higher across N. lapillus populations than equivalent S. umbilicalis populations. These findings agree with other studies comparing direct and indirect development as reproductive strategies across marine species (e.g., Sherman et al. 2008; Hoffman et al. 2011b), although reports of contrasting findings also exist (Ayre & Hughes 2000; Richards 2007). Putative outlier locus datasets identified more extensive genetic structure in both species than was identified by neutral genetic loci, suggesting a role for adaptive divergence despite ongoing gene flow. These results contribute to a growing body of literature showing that high gene flow and adaptive divergence can co-occur in marine organisms (e.g., Gleason & Burton 2016; Hoey & Pinsky 2018; Sandoval-Castillo et al. 2018; Xuereb et al. 2018; Selmoni et al. 2020). Seascape genomic modelling indicated a greater role of environment on genetic structure in N. lapillus than S. umbilicalis , with the latter being more attributable to spatial geographic patterns that coincide with the seasonal occurrence of oceanic fronts in the Irish Sea.
As hypothesized, patterns of genetic structure indicated greater connectivity across our study area in the broadcast spawning S. umbilicalis than the direct developer N. lapillus , regardless of the size of habitat gaps between sites. Although there were significant but low levels of differentiation (FST) between many sites and a marginally significant effect of IBD, there was a lack of distinct neutral genetic structure across the sampled area for S. umbilicalis . Larval connectivity (AEM vectors), geographic structure (dbMEM vectors), and environmental variables (PCs representing summer AT and SST, and wave exposure) best explained the variation in this dataset, with the largest proportion of explainable variation attributed to geographic structure. However, when partitioned, no partial RDA models were significant, supporting a general lack of explainable structure in the dataset. Taken together, the genetic homogeneity ofS. umbilicalis populations within our study area can at least partially be attributed to a combination of predicted larval dispersal over small and medium distances and multi-generational stepping-stone dispersal over larger distances. Surprisingly, predicted larval connectivity accounted for only a small proportion of variation in the neutral genetic dataset, perhaps because we did not account for interannual variability or climate change in larval connectivity matrices. Although Coscia et al. (2020) show interannual variability to be low over short time-scales (6 years), incorporating this variability may become more important over longer time-scales like those considered here (100 years). Further, what we considered to be large habitat gaps in our sampling design seemingly only constitute medium-scale dispersal distances for S. umbilicalis , given that genetic connectivity across these gaps was high.
Although we found much more genetic structure in N. lapillus thanS. umbilicalis , we found greater connectivity between populations separated by large habitat gaps than would be expected for a direct developing species (e.g., sites 2/3 and 10; Fig. 3). This finding supports the contribution of drifting/rafting (e.g., Colson & Hughes 2004), and stepping-stone dispersal (e.g., Crandall et al. 2012) to population connectivity in this species. Neutral genetic structure inN. lapillus indicated 5-11 genetic clusters of geographically proximate sites, and while there was a significant relationship between shortest marine coastal distance and FST, no significant effect of IBD was detected. Neutral genetic structure was best explained by two environmental PCs representing winter AT and summer AT and SST. The effects of environmental gradients (e.g., SST, AT) on putatively neutral SNP loci may be explained by either isolation by adaptation (IBA; Nosil et al. 2009), whereby strong ecological selection against immigrants results in adaptive population divergence that restricts gene flow and allows neutral loci to diverge via genetic drift (Thibert-Plante & Hendry 2010), or loose linkage to loci under selection (Gagnaire et al. 2015). Additionally, total explainable variation was low for this dataset, supporting the contribution of additional factors to describing neutral genetic structure in this species. Future efforts would benefit from incorporating simulations of putative long-distance dispersal of N. lapillus (e.g., rafting).
Interestingly, sites within habitat gaps were more likely to be colonized by N. lapillus , which we often found in small patches or sub-optimal habitat that did not support populations of S. umbilicalis . This may suggest that N. lapillus is more of a habitat generalist than S. umbilicalis , or, alternatively, that direct developers will have a greater likelihood of establishing new populations since founders can be a fertilized female or a drifting egg mass, allowing multiple offspring to hatch within the same area (Johannesson 1988). This tactic facilitates rapid population increase as encounter rates between individuals will be high in species with low mobility and high rates of self-recruitment. In contrast, while S. umbilicalis was found at two sites within the north Wales–Scotland gap, we found only one individual at one of these sites (St. Bees Head), and a very small population in the other (Selker Bay). This suggests that while larval dispersing species may reach distant sites more quickly and frequently than direct developers, they are likely to be spread over a much broader area during their planktonic larval feeding stage (Johannesson 1988). This dispersal tactic can lead to a low density of individuals at sites beyond a critical distance, and thus low population sizes and reduced encounter rates for subsequent generations of reproduction.
In contrast to a lack of neutral genetic structure in S. umbilicalis , outlier loci distinguished four genetic clusters indicating that selection can still be a major driver of spatial genomic structure, even in the face of extensive gene flow. These results are not surprising for a broadcast spawning species, for which large populations sizes are expected to reduce the effects of genetic drift and thus increase the probability that population differentiation results from localized natural selection (Nielsen et al. 2009; Gagnaire et al. 2015). The distribution of adaptive genetic structure in S. umbilicalis corresponds to clusters separated by the habitat gaps we investigated. This structure was significantly explained by both geographic and larval dispersal variables, but none of the environmental variables investigated. However, larval dispersal vectors were not significant in partial RDA models. Rather, most of the variation was attributed to geographic vectors representing large-scale spatial patterns differentiating sites in the northern Irish Sea/North Channel from all others (Fig. 6, dbMEM_2), and Irish from British (particularly south Wales) sites (Fig. 6, dbMEM_3). Interestingly, patterns in the spatial predictor dbMEM_2 closely resemble the “Forbes Line,” which demarcates the general northern limit of southern species within the UK and Ireland (Forbes 1858), and in the summer coincides with tidal fronts in the Irish Sea (Simpson & Hunter 1974). Where they occur, fronts have been shown to affect the dispersal (e.g., Ayata et al. 2010; Firth et al. 2021), and survival (e.g., Gaylord & Gaines 2000) of species by establishing spatially stratified environmental conditions in temperature and/or salinity (Pingree et al. 1974; Pineda 1994). The extent to which these fronts may drive adaptation in S. umbilicalis at its northern range edge is unknown, but seasonal fronts that form in the Irish Sea during summer (Simpson et al. 2009) may be especially influential, as they coincide with the timing of spawning and dispersal. Alternatively, these results may suggest that our models did not include other important environmental variables driving selection inS. umbilicalis or that selective pressures may vary spatio-temporally, resulting in patterns of chaotic genetic patchiness that do not easily correlate with environmental features.
The outlier locus dataset for N. lapillus suggested the potential for 7-8 adaptive genetic clusters that do not correspond completely with the structure observed in the neutral dataset. This putative adaptive structure was best described by the same environmental variables as for the neutral dataset (winter AT and summer AT and SST) but with the addition of wave exposure and a minor, non-significant contribution of one geographic vector. The significance of temperature to the adaptive structure of N. lapillus populations in our study area corroborates previous evidence supporting thermal-mediated selection in this species. Specifically, Chu et al. (2014) identified several fixed single nucleotide polymorphisms (SNPs) within heat stress-mediated genes between clades of N. lapillus in the northwestern Atlantic. Additionally, the significance of wave exposure corroborates many previous studies establishing relationships between exposure and adaptation of shell morphology in N. lapillus populations from Europe (Hughes & Taylor 1997; Guerra-Varela et al. 2009; Pascoal et al. 2012; Carro et al. 2019) and North America (Etter 1988, 1996).
The most supported driver of environmental adaptation in marine species is temperature (Liggins et al. 2019), and in the present study we provide further evidence to substantiate its important role. Temperature is as a strong stressor in the intertidal, where it underlies many important ecological (e.g., latitudinal distributions [Helmuth et al. 2006] and vertical zonation [Somero 2002]) and physiological (Tomanek & Helmuth 2002) processes. Intertidal species are influenced by both sea and air temperatures during high and low tides, respectively. Unfortunately, detailed and consistent meteorological records are generally unavailable for intertidal regions, where sea temperatures are likely influenced by air and ground temperature, substrate, aspect, and tidal range, among other factors. For this study we used average estimates of SST and AT over large, buffered areas (10 km), and thus acknowledge that we have likely not been able to characterize fine-scale environmental features of the local seascape that may substantially influence adaptation and fitness in our species. This may partially explain the lack of environmental associations detected in our S. umbilicalis outlier dataset, despite the known impacts of temperature on its dispersal and physiology.
S. umbilicalis spawns multiple times throughout the year in its range centre where sea temperatures are warmer, with shorter breeding periods towards its northern range edge (the current study area) where recruitment failure is associated with cooler temperatures (Bode et al. 1986; Kendall & Lewis 1986). Given the importance of temperature on larval dispersal and settlement in S. umbilicalis , fine-scale temperature data may aid in identifying environmental associations with outlier loci that are involved in reproductive processes.
Our inability to identify sequence matches for most of our outlier loci is not surprising, given the distant relationship of N. lapillusand S. umbilicalis to most model species, and the underrepresentation of genomic resources and sequencing efforts for marine invertebrates (Lopez et al. 2019). Of loci that we were able to match to NCBI sequences, few were within coding regions, an increasingly common result in SNP-mapping studies. Indeed, SNPs can be as far as two Mbp away from, and are not necessarily closest to the genes they affect (Brodie et al. 2016). Of the few loci that were located within coding regions, both species showed matches to the estrogen receptor (ER) gene. Although the function of the ER gene in molluscs is still largely uncharacterized, it has been shown to overexpress in N. lapillusin the presence of estrogenic chemicals in raw urban/industrial effluent pollutants, with concomitant increases in reproductive maturation (Castro et al. 2007). The other two outlier loci that mapped to coding regions in S. umbilicalis are less well characterized but may also be involved in reproduction (He et al. 2015; Robay et al. 2018).
Our results support the hypothesis that substantial heterogeneity of the seascape can support connectivity of the otherwise low mobility direct developing species, Nucella lapillus , and create adaptive divergence in an otherwise highly connected meta-population of the broadcast spawning species, Steromphala umbilicalis . Ultimately, characterizing spatial patterns in connectivity, estimating genetic diversity, and identifying locally adapted populations will lead to a better understanding of the resilience of marine species to changing environments and allow for improved management of fisheries and marine protected areas (MPAs; e.g., Miller & Ayre 2008; Sinclair-Waters et al. 2018). Sea temperatures around the UK and Ireland are predicted to catch up with global trends within the next decade, and associations between species abundance and sea temperatures observed over less than a decade indicate the sensitivity of many species to these changes (Mieszkowska et al. 2020). Thus, understanding the impacts of ocean warming on intertidal marine species, and their potential to adapt to these changes is an important goal.
Acknowledgements: This research was funded by the Ecostructure project (part-funded by the European Regional Development Fund (ERDF) through the Ireland-Wales Cooperation Programme 2014-2020). The authors would like to acknowledge all members of the Ecostructure project (http://www.ecostructureproject.eu/) who helped in the sample collection and processing for this work, and Data2Bio LLC (Ames, Iowa) who provided the sequencing services for this project. We are grateful to Niall McKeown, Peter Lawrence and Ally Evans for helpful discussions of the data and results. Data analysis was supported by the Supercomputing Wales programme, which is part-funded by the European Regional Development Fund (ERDF) via the Welsh government.