Fine‐scale seascape genomics of an exploited marine species, the common cockle Cerastoderma edule, using a multimodelling approach

Abstract Population dynamics of marine species that are sessile as adults are driven by oceanographic dispersal of larvae from spawning to nursery grounds. This is mediated by life‐history traits such as the timing and frequency of spawning, larval behaviour and duration, and settlement success. Here, we use 1725 single nucleotide polymorphisms (SNPs) to study the fine‐scale spatial genetic structure in the commercially important cockle species Cerastoderma edule and compare it to environmental variables and current‐mediated larval dispersal within a modelling framework. Hydrodynamic modelling employing the NEMO Atlantic Margin Model (AMM15) was used to simulate larval transport and estimate connectivity between populations during spawning months (April–September), factoring in larval duration and interannual variability of ocean currents. Results at neutral loci reveal the existence of three separate genetic clusters (mean F ST = 0.021) within a relatively fine spatial scale in the north‐west Atlantic. Environmental association analysis indicates that oceanographic currents and geographic proximity explain over 20% of the variance observed at neutral loci, while genetic variance (71%) at outlier loci was explained by sea surface temperature extremes. These results fill an important knowledge gap in the management of a commercially important and overexploited species, bringing us closer to understanding the role of larval dispersal in connecting populations at a fine geographic scale.

population differentiation within the marine realm (Benestan et al., 2015;Maroso et al., 2018). Nevertheless, there has been a shift in focus from neutral variation to adaptive genetic differentiation when studying marine organisms (Nielsen, Hemmer-Hansen, Larsen, & Bekkevold, 2009). Markers under selection display higher divergent frequencies and can facilitate the identification of population structure and individual assignment (Araneda, Larraín, Hecht, & Narum, 2016;Nielsen et al., 2012;Woodings et al., 2018). They can also indicate the potential for resilience to environmental change (Razgour et al., 2019). Connectivity between marine populations is central to their health and resilience to external pressures such as parasites and pathogens (Rowley et al., 2014), pollution, human exploitation and climate change over ecological and evolutionary timescales (Burgess, Bowring, & Shen, 2014;Cowen & Sponaugle, 2009;Gimenez, 2019). For these reasons, it is vital to distinguish between neutral variation and adaptive divergence when attempting to understand what drives the observed population structure in marine organisms.
A seascape genomics approach is particularly valuable in this context (Grummer et al., 2019;Selkoe et al., 2016). Genetic and genomic data can be used in conjunction with environmental variables such as sea water temperature, salinity, water depth, irradiance, turbidity and sediment type (Viricel & Rosel, 2014) against neutral and adaptive genetic differentiation (Coscia, Robins, Porter, Malham, & Ironside, 2012;Young et al., 2015). This can provide new insights for interpreting genetic patchiness in relation to specific environmental features (Benestan et al., 2016;Bernatchez et al., 2019;Truelove et al., 2017).
Seascape genomics has proven to be particularly useful when considering exploited species, as it has the potential to inform management and aid sustainable exploitation by enabling management units to be defined (Bernatchez et al., 2017;Teacher, André, Jonsson, & Merilä, 2013). Despite intense exploitation of many marine species, their management rarely takes into account genomic information (Bernatchez et al., 2017;ICES, 2018) and exploited aquatic invertebrates (shellfish) in particular receive little attention from policymakers and stakeholders in comparison with fish (Elliott & Holden, 2017). In the Irish Sea, shellfish represent 80% of the total landings per year, with an economic value of £46.6 million (Elliott & Holden, 2017). Among these, the common cockle Cerastoderma edule fisheries are some of the most valuable fisheries for both Ireland and Britain (Dare, Bell, Walker, & Bannister, 2004;Hervas, Tully, Hickey, O'Keeffe, & Kelly, 2008) and are of high socio-economic importance, valued at £3.3 million for Wales alone (Elliott & Holden, 2017).
The common cockle has both ecological and commercial importance, providing an important food source for wading birds in addition to employment for coastal communities (Flach & de Bruin, 1994;Hickin, 2013). Cerastoderma edule occurs in intertidal soft sediment regions of the eastern Atlantic, from Norway to Senegal. It can live for up to ten years and is characterized by high fecundity and high dispersal potential due to a pelagic larval phase which lasts for approximately 3-5 weeks following spawning from May to August (Malham, Hutchinson, & Longshaw, 2012). In recent decades, cockle stocks have declined across Europe, with production falling from 108,000 tons in 1987 to less than 25,000 tons in 2008 (Martínez, Mendez, Insua, Arias-Pérez, & Freire, 2013). Cockle declines have been attributed to different factors in different locations, such as overharvesting (Wolff, 2005) and parasitic infections (Longshaw & Malham, 2015;Thieltges, 2006). In the UK, recurrent mass mortalities have occurred at several long-established cockle fisheries, resulting in significant economic losses (Woolmer, 2013). These mortalities have not been attributed to any single environmental factor, and interactions between multiple stress factors are suspected (Callaway, Burdon, Deasey, Mazik, & Elliott, 2013;Malham et al., 2012). Sustainable management of the common cockle is hindered by poor understanding of their population connectivity. Analysis of microsatellite and mitochondrial DNA markers suggests weak barriers to gene flow between C. edule populations along the North European coast Martínez, Freire, Arias-Pérez, Mendez, & Insua, 2015). However, these markers lack sufficient resolution to investigate connectivity at the finer scales relevant to fisheries management (Bernatchez et al., 2017).
Given the major logistical challenges of directly quantifying larval connectivity, efforts have focused on simulating ocean hydrodynamics to estimate larval dispersal (Cowen, Gawarkiewicz, Pineda, Thorrold, & Werner, 2007;Paris, Cowen, Claro, & Lindeman, 2005;Robins, Neill, Giménez, Jenkins, & Malham, 2013). This approach identifies well-connected population groups as well as weakly connected, partially connected or isolated populations. These simulations highlight the importance of local and mesoscale hydrodynamics interacting with species-specific larval behaviours in driving population persistence (Bode et al., 2019;Botsford et al., 2009;North et al., 2008;Robins et al., 2013) and recovery from stock decline (Gimenez, 2019). They highlight how the capacity of a population to recover from mass mortalities is contingent on the scale of disturbance relative to the scale of connectivity (Masier & Bonte, 2019). Usual circulation patterns, and hence connectivity, can be modulated by severe wind and wave conditions. Previous larval dispersal studies have predicted that given atypical meteorological conditions during spawning events, new connectivity routes can be established (e.g. by reversing the Celtic Sea front circulation (Hartnett, Berry, Tully, & Dabrowski, 2007) or by large-distance displacements from expected routes affecting sea turtles (Monzón-Argüello et al., 2012)).
In the present study, a seascape genomics approach using single nucleotide polymorphisms (SNPs) was employed, for the first time, to resolve patterns of population structure of the common cockle between estuaries within a commercially active area (the Irish and Celtic seas), with the goal of identifying management units. For this purpose, we first tested for neutral population structure and then assessed the role of current-mediated larval dispersal in shaping it.
We then investigated the relationship between environmental factors and adaptive divergence, after identifying outlier markers, with a particular focus on abiotic factors such as water temperatures (surface) and ocean stratification during the spawning season. To do this, larval transport between sites was estimated and connectivity matrices derived from oceanographic modelling, accounting for interannual variability due to biophysical parameters, like spawning time and larval duration.

| RAD sequencing and bioinformatic analysis
Reduced representation libraries (Baird et al., 2008) were constructed using the restriction enzyme PstI (New England Biolabs) for restriction site-associated DNA sequencing (RADseq). RAD libraries (each consisting of 24 uniquely barcoded individuals) were produced according to Etter, Preston, Bassham, Cresko, and Johnson (2011).
Each library was quantified using real-time PCR and single-end (100bp target)-sequenced in-house on an Illumina HiSeq 2000.
Initial bioinformatic analysis, including quality control, demultiplexing and identification of polymorphisms, was performed by Floragenex (www.flora genex.com), using Samtools 0.2 (Li et al., 2009) and custom scripts, retaining one SNP per tag, with a minimum coverage depth of six for each allele, and genotyped in at least 70% of the individuals in the overall data set. These settings produced an initial data set of 191 individuals and 4,271 single nucleotide polymorphisms (SNPs), which was then further filtered by the authors using the packages poppr 2.8.1 (Kamvar, Brooks, & Grünwald, 2015;Kamvar, Tabima, & Grünwald, 2014) and adegenet (Jombart, 2008;Jombart & Ahmed, 2011) in R (R Core Team, 2019). Markers with data missing in more than 25% of individuals were discarded from the data set, as well as loci with F IS equal to 1, −1 or NA. Three MAF (minimum allele frequency) filters of 0.05, 0.025 and 0.01 were then applied, generating three data sets. No significant effect of MAF filter upon heterozygosity, global F ST and population structure was detected, and so a MAF filter of 0.01 across all sites was applied (Xuereb et al.., 2018). This allowed us to maximize the number of markers and hence the information available at fine spatial scale, while reducing the bias that might be introduced by retaining low-frequency SNPs (Roesti, Salzburger, & Berner, 2012). All the downstream analyses were thus performed on the MAF = 0.01 data set. Markers in linkage disequilibrium (LD) were identified and removed using the function pair.ia of the poppr package (r 2 > .7). Finally, SNPs that were found to deviate from Hardy-Weinberg equilibrium (HWE; p = .01) in four out of seven populations were removed (Wyngaarden et al., 2018).
The final, filtered data set included 138 individuals from the seven locations, and 1,725 SNPs (Table S1).

| Neutrality tests and population structure
Genomic markers were tested for neutrality using two complementary methods: BayeScan 2.1 (Foll & Gaggiotti, 2008) and the R package pcadapt (Luu, Bazin, & Blum, 2017). In order to minimize the detection of false positives, we considered as outliers only those SNPs that were selected by both methods . BayeScan is an F ST -based, Bayesian method (Beaumont & Balding, 2004), while pcadapt is based on a principal component analysis (PCA) of individual genotypes and is known to perform particularly well in the presence of weak structure, admixture or range expansions (Luu et al. 2017).
BayeScan was run using default settings (5,000 iterations, 10 thinning intervals, 20 pilot runs-5,000 iterations each, 10 prior odds). For pcadapt, we chose the most appropriate number of clusters K in the scree plot, which displays in decreasing order the percentage of variance explained by each principal component (PC). Q-values were used to account for false discovery rate, and SNPs were considered as significant outliers at alpha values ≤.05. Genetic diversity was estimated on three data sets: overall, neutral and outliers, in order to disentangle the role of demographic processes versus selection. Heterozygosity (expected and observed) was estimated with the R package hierfstat (Goudet, 2005;Goudet & Jombart, 2015), while population pairwise F ST (Weir & Cockerham, 1984) and relative 95% confidence interval (1,000 bootstraps) were estimated with the R package assigner (Gosselin, 2019).
Individual-based population structure was assessed using two approaches: discriminant analysis of principal components (DAPC) as implemented in adegenet (Jombart, 2008;Jombart & Ahmed, 2011) and the function sNMF of the LEA package (Frichot & Francois, 2015). While DAPC uses a priori spatial information (sampling locations), sNMF does not make any assumptions about sampled and ancestral populations. The number of K clusters in the data sets was chosen using the Bayesian information criterion (Schwarz, 1978) in DAPC. In sNMF, the entropy criterion provided a basis for choosing the number of ancestral populations that best explain the genotypic data (Alexander & Lange, 2011;Frichot, Mathieu, Trouillon, Bouchard, & Francois, 2014).
The power of the neutral and outlier data sets to discriminate and assign individuals was determined with a genotype accumulation curve ( Figure S2), calculated with the function genotype_curve of the poppr package.

| Hydrodynamic modelling
Simulations of 3D flow fields were used to drive a particle tracking model (PTM) developed to simulate potential larval transport and F I G U R E 1 (a) Sampling locations; (b) DAPC using the neutral data set; (c) DAPC using the 14 outliers; below, the barplots generated by sNMF from the neutral (d, f) and outlier (e, g) data sets, for both K = 2 and K = 3

| Particle tracking model (PTM)
The particle tracking model simulates the lagrangian movement of individual particles in space and time, based on oceanographic dispersion and individual particle behaviour. In our case, particle behaviour was represented by varying the pelagic larval duration to reflect the observations in the field . However, given the paucity of data about C. edule larval behaviour, the effects of vertical swimming were not incorporated into the PTM. The PTM was programmed in Python.
A previous sensitivity study by Robins et al. (2013)  Values that are significant (95% confidence interval) are in bold.
TA B L E 1 Genetic diversity indices for the two data sets.
F I G U R E 2 (Left) Simulated mean depth averaged ocean current strength in the Irish and Celtic seas for April (top) and July (bottom) 2014. The magnitude of the currents is shaded in red, and the direction of the currents is indicated with the grey arrows. (Right) Simulated mean surface to bottom temperature difference for April (top) and July (bottom) 2014. This is a measure of the degree of stratification of the water column with blue colours resembling a more mixed water column and yellows a more stratified water column were assumed to have settled there. The 10 km threshold was based on the average tidal excursion for the Irish Sea (Robins et al., 2013), although a range of thresholds are discussed in the Section 3.

| Spatial eigenfunction analysis (SEA) and environmental association (EA) analysis
To test for association between the oceanographic environment and the cockle genetic structure, a spatial eigenfunction analy- parsimonious RDAs were carried out using the variables selected (Borcard, Gillet, & Legendre, 2011).

| Genetic diversity and population structure
After filtering for missing data, F IS , minimum allele frequency,   Genotype accumulation curves ( Figure S2) for both data setsneutral and outlier-show that a plateau is reached and variance decreased, indicating that the data sets have sufficient power to F I G U R E 3 (a-g) Probability density distribution maps for 2014 showing simulated dispersal probability from release Sites 1-7 (black squares). Each panel shows dispersal probability for 11,520,000 particles (12,000 each month × 6 months × 10 settlement days). (h) Seasonally-averaged connectivity networks: the thickness of the pathways in corresponds to the average modelled connectivity. The size of the red circles corresponds to self-recruitment, and the dashed green circle show the settlement radius used to estimate connectivity distinguish between individuals (Arnaud-Haond, Duarte, Alberto, & Serrão, 2007).

| Seasonal circulation and larval dispersal
The simulated circulation structure within the Celtic and Irish seas has, in our study, contributed markedly to distinct patterns of larval transport and population connectivity that would not be anticipated from geographic location alone. Moreover, spatial variability in these patterns was simulated over seasonal timescales. The strong tides that characterize this region and have the potential to advect larvae tens of kilometres away but are oscillatory, resulting in minimal net larval dispersal (Robins et al., 2013;Robins et al. 2015). Thus, larval transport is controlled by density-driven currents (or lack of) along frontal boundaries that develop due to thermal stratification during summer months (Horsburgh, Hill, & Brown, 1998;Simpson & Hunter, 1974). Although much weaker than tidal currents, these density-driven currents are persistent over time and act as key pathways for dispersing larvae. For example, the Celtic Sea front may facilitate connectivity of cockle populations from South Wales across to Ireland , but also restrict transport across the front between the Celtic and Irish seas.  (Figure 3 and Figure S5). This is shown by comparing the simulated particle dispersal from Dee (Site 1, enclosed bay) with  (Figure 3a and Figure S5). In contrast, particles from Red Wharf Bay were exposed to stronger currents and the site experienced lower self-recruitment (13.6 ± 3.5%) with "hot spots" of higher density particles (up to 0.5%) advected ~100 km northwards ( Figure 3b).
Our simulations suggest a high degree of isolation at Dyfi (Site 3; Figure 3c and Figure S5), with the majority of particles retained in close proximity. At the same time, due to the weak flows, simulated particles from other sites never reached this site. Particles from Gann (Site 4) and Burry (Site 5) were capable of dispersing readily between one another (up to 1% connectivity) and also westwards to the Irish populations, particularly from the more exposed Gann Estuary (up to 0.5% connectivity) (Figure 3d,e and Figure S5). Connectivity with Dee was possible but unlikely (<0.01%
Partial RDAs on neutral data using these three variables were not significant.
These SST values were the lowest daily mean SST in each month ( Figure S1). The global parsimonious RDA (Figure 4a) was overall nonsignificant when including the three selected factors at once (p = .11), although it was significant when including two at a time: geographic distance (dbMEM-1) + July mean SBTD (p = .008), and simulated connectivity (AEM2) + geographic distance (dbMEM-1) (p = .022). Partial RDAs were not significant. The parsimonious RDA run on outliers was globally highly significant (p = .002; Figure 4b).
For the neutral RDA, the first axis explained almost 48% of the total variance and was mainly driven by geographic distance

| D ISCUSS I ON
This is the first study to use genomic markers and biophysical lar- Overall, the genetic results fit well with the predictions of the larval transport model, providing a level of empirical validation for both the simulated hydrodynamics and connectivity by larval dispersal. Environmental association analysis revealed that neutral genetic structure was strongly linked to geographic distance between sites and to the strength and direction of the ocean currents acting as corridors for larval dispersal, whereas colder periods (cold SSTs) were identified as potential drivers of adaptive divergence.

| Connectivity, fine-scale population structure and adaptive divergence
Neutral genetic diversity was low across the cockle populations within the study area, and low compared with previous studies on marine bivalve population genomics (e.g. Bernatchez et al.,

2019; Lehnert et al., 2019). Populations of cockles in the Irish
Sea have been under pressure for at least two decades, with mass mortality events and declines due to overexploitation leading to strict management of most beds across the UK (Woolmer, 2013). In addition, variance in reproductive success is known to occur in bivalves (Hedgecock & Pudovkin, 2011). These events could be responsible for the loss of genetic diversity, as already observed in several marine organisms (Pinsky & Palumbi, 2014).
On the other hand, heterozygote deficiency (as indicated by positive F IS in all sampling locations; Table 1) is well known to occur in F I G U R E 4 Redundancy analysis (parsimonious RDA) performed for the (a) neutral and (b) outlier data sets. Each circle is a sampling location, and each arrow is an environmental variable that significantly drives the observed population structure marine bivalves (Gaffney, 1994), and had already been detected in these cockle populations using microsatellites .

Furthermore, cockles have been shown to undergo boom and bust
years (Morgan, O'Riordan, & Culloty, 2013) with the dispersal of cockle larvae and recruitment altering between years (Miller, Versace, Matthews, Montgomery, & Bowie, 2013;Morgan et al., 2013). Another explanation of positive F IS could be the Wahlund effect, caused by the presence of substructure within populations.
This was not considered to be the case in this data set, given that no subpopulation pattern was found in the clustering analysis (Manangwa et al., 2019).
Given the fine spatial scale and the reproductive biology of the study organism (broadcast spawner with a long pelagic larval phase), a lack of or weak population structure was expected.
Geographic proximity certainly favours gene flow, but oceanographic currents aiding larval transport are also major drivers of population structure (Barbut et al., 2019). For example, Red Wharf Bay and the Dee Estuary appear to be genetically homogeneous, due to high levels of gene flow, but these populations are very distinct from those further south (150-350 km away). This concurs with the model's prediction that larvae from the north coast of Wales will disperse northwards. The high levels of gene flow detected between the Gann Estuary on Wales' south-west coast and sites on the southeast coast of Ireland also corresponds with the model's prediction of high levels of westward dispersal along the Celtic Sea front .
Of particular interest is the genetic make-up of the population in the Burry Inlet. Here, cockles have experienced recurrent mass mortality events for 15 years (https://marin escie nce.blog. gov.uk/2015/08/14/unusu al-cockle-morta lities-burry-inlet/). The relatively high level of genetic differentiation between the neighbouring Burry Inlet and Gann Estuary populations indicates that gene flow between these populations is low, seemingly contradicting the larval transport model's prediction of high connectivity. This result suggests that the Burry Inlet population may able to persist through self-recruitment, rather than forming a sink population depending upon immigration from healthier populations elsewhere, and could explain the low genetic diversity and high levels of inbreeding detected there. However, the model showed that larval dispersal is possible from the Gann to Burry estuaries, although we assumed a relatively large settlement zone that extended beyond the mouth of the Burry Inlet ( Figure S9).
Furthermore, it must be acknowledged that the Gann Estuary cockle population is far smaller than that of the Burry Inlet so larvae dispersing from the Gann into the Burry may be swamped by self-recruitment or might not survive due to local selection against them (i.e. density dependence; Ford, Shima, & Swearer, 2016).
Additionally, the southern coast of Wales contains other large cockle populations, such as the Three Rivers fishery midway between the Gann and Burry. These populations were not sampled/ modelled in this study but may provide a nearer and greater source of larvae for the Burry Inlet than does the Gann Estuary.
The spatial environmental association analysis identified several environmental factors that appeared associated with the population structure observed at neutral and outlier genetic markers. This is a strong statistical approach, which has already been successfully employed to study the influence of the environment on the genetic structure of commercially important bivalves, such as the eastern oyster (Crassostrea virginica; Bernatchez et al., 2019) and Atlantic deep-sea scallop (Placopecten magellanicus; Lehnert et al., 2019).
In C. edule, neutral genetic structure is strongly dependent on geographic distance between sites (dbMEM-1), indicating that isolation by distance plays an important role in shaping the observed genetic structure in this species, despite its long pelagic larval duration.
Nevertheless, it is the interplay between isolation by distance, extremes in temperature-driven currents (July SBTD) and simulated connectivity (AEM2) that shapes neutral population structure.
The summer stratification which strengthens the Celtic Sea front current, directed from south Wales to Ireland (Simpson & Hunter, 1974), Furthermore, our connectivity matrices (Figure 3 and Figure   S5) were not able to take into account the population abundance of cockles and hence the number of larvae produces at each site.
Our modelling results suggest that larval dispersal is stable on interannual timescales (2008)(2009)(2010)(2011)(2012)(2013)(2014). However, a degree of interannual variability was detected along exposed coasts and in re- Considering the 14 outlier loci, the minimum sea surface temperatures recorded in April and in September 2014 explained most of the genetic variance observed within this data set. In particular, the Burry Inlet was strongly associated with the SSTs in September, which was warmer compared with SSTs recorded at other locations for the same time ( Figure S1). If cockles in the Burry Inlet are indeed adapted to warmer sea surface temperatures than populations at sites further west or north (and especially compared with SSTs at Gann), it could be speculated that this may explain the maintenance of genetic differentiation between the Burry Inlet and the Gann Estuary despite their spatial proximity and the potential for high connectivity predicted by the model (Figure 3 and Figure S5). Larvae dispersing from Gann to the Burry may not be able to survive postrecruitment given selective pressure against them by higher temperatures. The data collected in this study is not suitable to test this hypothesis, but future investigations of local adaptation of bivalves in the Irish Sea should take these findings into account. Because of the lack of a reference genome and because the sampling locations are not representative of an environmental gradient, our data cannot be used to test this hypothesis of adaptation.

| Implications for management
The results from our modelling are likely to be relevant to fishery management in terms of the potential seasonal and interannual variability in larval supply to cockle grounds. This study demonstrates the existence of three distinct units of cockles using both genomic tools and larval dispersal modelling. As with other recent studies (Lal, Southgate, Jerry, Bosserelle, & Zenger, 2017) these findings have important implications for fishery management Miller et al., 2013) and how fisheries management can be reconciled with conservation and other activities.
Given the incidence of recurrent mass mortality events at the Burry Inlet (Callaway et al., 2013), the genetic isolation of this cockle fishery implied by the results of this study should be investigated further. This could be achieved by expanding the sampling coverage of the Burry Inlet to multiple sites and years and assessing its connectivity to other nearby cockle beds that have not been included in this study. Future investigations should be aimed at clarifying the role of adaptive divergence into the fine-scale population dynamics of the common cockle in this area, to improve management while also assessing the role played by diseases and infections. The results from this study highlight the importance of the use of genomic and hydrodynamic data in assessing population structure and connectivity in an exploited and commercially important marine species and may aid in current and long-term management regimes of this species (Bernatchez et al., 2017;Lal et al., 2017).

ACK N OWLED G EM ENTS
The genomic data were generated thanks to an internal grant awarded to IC by Aberystwyth University. The authors also wish to acknowledge the support of the Interreg Ireland-Wales Programme ISPP

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data sets used in this study are available from the Dryad Digital Repository (https://doi.org/10.5061/dryad.ht76h drbr).