- Research article
- Open Access
Seascape genomics reveals population isolation in the reef-building honeycomb worm, Sabellaria alveolata (L.)
BMC Evolutionary Biology volume 20, Article number: 100 (2020)
Under the threat of climate change populations can disperse, acclimatise or evolve in order to avoid fitness loss. In light of this, it is important to understand neutral gene flow patterns as a measure of dispersal potential, but also adaptive genetic variation as a measure of evolutionary potential. In order to assess genetic variation and how this relates to environment in the honeycomb worm (Sabellaria alveolata (L.)), a reef-building polychaete that supports high biodiversity, we carried out RAD sequencing using individuals from along its complete latitudinal range. Patterns of neutral population genetic structure were compared to larval dispersal as predicted by ocean circulation modelling, and outlier analyses and genotype-environment association tests were used to attempt to identify loci under selection in relation to local temperature data.
We genotyped 482 filtered SNPs, from 68 individuals across nine sites, 27 of which were identified as outliers using BAYESCAN and ARLEQUIN. All outlier loci were potentially under balancing selection, despite previous evidence of local adaptation in the system. Limited gene flow was observed among reef-sites (FST = 0.28 ± 0.10), in line with the low dispersal potential identified by the larval dispersal models. The North Atlantic reef emerged as a distinct population and this was linked to high local larval retention and the effect of the North Atlantic Current on dispersal.
As an isolated population, with limited potential for natural genetic or demographic augmentation from other reefs, the North Atlantic site warrants conservation attention in order to preserve not only this species, but above all the crucial functional ecological roles that are associated with their bioconstructions. Our study highlights the utility of using seascape genomics to identify populations of conservation concern.
Conservation strategies in the context of climate change rest on the premise that populations can show three responses to avoid fitness loss: evade, evolve or acclimatise . Characterisation of genetic diversity can shed light on these responses , by indicating dispersal and gene flow using neutral genetic markers (evasion), or adaptation using adaptive genetic variation, which inform us about the potential for an organism to evolve in altered conditions [3,4,5]. Separating neutral and adaptive causes of phenotypic differentiation, or lack thereof in a heterogeneous environment (sensu countergradient variation ;), can require logistically difficult experiments . However, advances in genomic technologies and analyses have opened up new avenues for exploring genome-wide (neutral) versus locus-specific (adaptive) genetic variation between populations, allowing us to compare patterns of adaptive genetic diversity with environmental variation to infer agents of selection [8, 9]. These methods have already been successful in contrasting neutral and non-neutral genetic variation in a range of taxa including fish [5, 8, 10, 11], insects , and marine invertebrates [12, 13] and represent an important step in interpreting genetic diversity in relation to environment.
Despite genomic advances allowing investigations of genetic diversity in non-model species, which is particularly important for species of conservation concern, understanding of gene flow and local adaptation in marine environments had, until recently, lagged behind that for terrestrial species [14,15,16]. Recent advances in seascape genomics have facilitated consideration of the effect of landscape on neutral and adaptive genetic variation in marine environments, taking into account ocean circulation patterns and environmental variation to explain observed population structuring [8, 17]. For instance, Benestan et al. (2016)  found that ocean currents as well as geographic distance were key to explaining observed patterns of neutral genetic variation in the American lobster (Homarus americanus H. Milne Edwards, 1837), whereas sea surface temperature acted as an agent for natural selection. Ocean circulation modelling for understanding dispersing larval stages has emerged as an important tool to explain neutral population structure where isolation-by-distance measures alone are uninformative [18,19,20].
Intertidal communities have been proposed as particularly suitable for studying the impacts of climate change, as many species already exist at the limit of their thermal tolerance ranges [21, 22]. Ecosystem engineers provide and modify habitat for other species and are thus disproportionately important to diversity, ecological functioning, and conservation in a changing climate [1, 14, 23,24,25]. Therefore, understanding how littoral ecosystem engineers have adapted to the local environment they experience is key to understanding not only their future survival in a changing climate, but also the potential survival of intertidal ecosystems as a whole.
The honeycomb worm (Sabellaria alveolata (Linneaus, 1767)), is an ecosystem engineer that constructs biogenic reefs along temperate coastlines of North Africa and Europe . They are gonochoric sexually reproducing broadcast spawners. Asexual reproduction is not known to occur in the species, in line with the observation that clonality occurs in less than 1% of polychaetes . Larvae are dispersive but preferentially settle on pre-existing reef or remains, although they do not discriminate based on reef or origin [28, 29]. Each individual constructs a tube from sand, which is glued together with biomineralised cement [29, 30]. Individuals are gregarious and tubes aggregate to form large bioconstructions or reefs, which can stretch over several kilometres [31,32,33]. These reefs provide a large diversity of microhabitats for other species and alter wave action, water circulation, erosion and sedimentation patterns , and phytoplankton concentration  and thus modify the availability of resources, leading to a high associated biodiversity compared to surrounding areas [36, 37]. Biogenic reefs, including those of Sabellaria, are listed under the EC Habitats Directive (92/43/EEC) Annex 1 and are thus of conservation interest . Despite this, little is known about neutral genetic structure (but see [39, 40]) and to our knowledge nothing is known about adaptive genetic structure of S. alveolata. Neutral population structuring is expected because (i) adults are sessile, and (ii) during the dispersive planktotrophic larval stage of this species, individuals have been shown to demonstrate positive taxis towards the cement of sand grains and bioclasts of S. alveolata tubes  leading to site philopatry . However, modelling of larval dispersal is needed to fully understanding how the environment influences gene flow within marine systems .
We have previously demonstrated that S. alveolata show local adaptation to temperature in terms of their membrane lipid composition  but the genomic basis of this adaptation is not known. Therefore, as their distributional range extends along a latitudinal thermal gradient from Scotland to Morocco, S. alveolata provide an excellent opportunity to investigate if/how adaptive population structure is related to temperature, and assess the role of neutral processes in gene flow.
In this study, we explore the extent of neutral and adaptive genetic variation among reef sites of the honeycomb worm, S. alveolata, and the relationship with the environment. We test: 1) if gene flow is restricted among S. alveolata reef sites and whether population structure can be predicted by ocean circulation modelling; and 2) whether S. alveolata show evidence of local adaptation and if this relates to local thermal conditions.
SNP detection and summary statistics
Due to high variability in sequence coverage between samples (37,000 - 1,400,000 reads per sample; Additional file 1), only 68 samples across the nine sites were retained for analysis (Table 1). Of these, 50 samples were from sites where temperature data was also available (Table 1; Fig. 1). In total, 439,598 SNPs were identified across 273,397 RADtags with an average coverage depth of 14.85x per sample. Due to low read coverage overall, only 506 variable SNPs were identified, 482 of which were from independent RADtags and could be used for population genetic analyses. Marker independence across sites was confirmed with no loci flagged as showing linkage disequilibrium. Although some loci showed significant deviation from Hardy-Weinberg, the loci involved differed across sites and were thus included in further analyses .
All sites showed reasonable numbers of alleles per locus, nucleotide diversity and expected heterozygosity (N alleles = 2.44 ± 0.13, θ = 0.26 ± 0.05, He = 0.52 ± 0.04). Observed heterozygosity was low across all sites (Ho = 0.11 ± 0.02) but was significantly lower than He only in the Balearic Sea site, which had a low sample size (Table 1). Most molecular variance was found within rather than between sites (within = 71.42%, between = 28.58%, P < 0.01) and within rather than between the genetic clusters identified by STRUCTURE (see below) (within = 84.08%, between = 15.92%, P < 0.01).
From the total number of SNPs, ARLEQUIN identified 191 as outliers (39.63%) and BAYESCAN identified 28 as outliers (Fig. 2; 5.81%). Twenty seven of the 28 outlier SNPs identified by BAYESCAN as outliers were also identified by ARLEQUIN. Therefore, the 27 outlier SNPs (5.60%) were removed from the data set and 455 (94.40%) were used in the neutral genetic analyses.
Pairwise FST values for neutral loci had an average of 0.28 ± 0.10 across sites, with 25 out of the 36 comparisons showing significance following Bonferroni correction, strongly demonstrating a high degree of population substructure within the system as a whole (Table 2). In particular, the North Atlantic, Tyrrhenian Sea and Bay of Biscay sites were genetically differentiated from all other sites (Average FST = 0.42 ± 0.10; 0.32 ± 0.12 and 0.25 ± 0.05, respectively; Table 2), suggesting very low gene flow to and from these sites. The Discriminant Analysis of Principal Components (DAPC) results also showed the Tyrrhenian Sea to be highly differentiated from the other sites (Fig. 3a). When the Tyrrhenian Sea was excluded from the analysis, the English Channel population emerged as being differentiated from the remaining populations (Fig. 3b). Using the Bayesian cluster analysis in STRUCTURE, ΔK showed a clear peak at K = 2, suggesting two genetic clusters within the nine sites studied (Fig. 4a). Likelihood probability profiles also showed that K = 2 had the highest likelihood with low variance. The barplot in STRUCTURE for K = 2 identified that individuals from the North Atlantic site formed a distinct population from all the other sites combined (Fig. 4b). This result did not vary regardless of whether sampling location was included as a prior and there was no evidence of hierarchical clustering (K = 1 when North Atlantic individuals were removed). FST based isolation (Fig. 5a) by straight-line geographic distance (Table 2) and shortest ocean distance (Table 2; Fig. 5b) were not significant (r = 0.14, P = 0.27 and r = 0.08, P = 0.27, respectively).
Maps of larval densities from dispersal simulations show the range of potential dispersal pathways for a single generation released from each sample location: from the most conservative (1 week of 6 hourly releases, 5 weeks maximum planktonic larval duration (PLD), Fig. 6a), to the most optimistic (daily releases all year round, 12 weeks PLD, Fig. 6b). The most optimistic dispersal scenarios are suggestive of isolation among sites (Additional file 1). Model simulations suggest isolation among the North Atlantic, northern Irish Sea, and southern Irish Sea sites (Additional file 1), which is reinforced by the FST distances (Table 2). Under the conservative scenario estimates, no larvae were exchanged among any sites in a single generation (Fig. 6a). Under the most optimistic scenario, high densities of larvae were retained close to natal release sites (average within 133 km, maximum 230 km) and very few larvae exchanged among release locations (Additional file 1). Potential connectivity was only identified between the northern Irish Sea and southern Irish Sea sites (as per ) but < 0.1% of the larval population made this transition. Dispersal distances away from source were shortest from the North Atlantic release site on the west coast of Ireland; this site exhibited the highest proportion of larval retention (79% self-recruited) and very few larvae reached nearby reefs, and none reached other release sites (Additional file 1). Dispersal resistance, as identified using ocean circulation modelling, did not show significant correlation with genetic distance, in terms of FST, for any of the PLD scenarios (5 weeks: r = − 0.02, P = 0.59; 8 weeks: r = − 0.01, P = 0.59; 10 weeks: r = − 0.006, P = 0.57; 12 weeks: r = 0.00, P = 0.57).
Loci under selection
Both ARLEQUIN and BAYESCAN analyses suggested that the 27 loci identified as outliers were under balancing selection (Fig. 2). Temperature data could not be collected from three of the study sites due to loss of equipment during the data collection period (Balearic Sea, Tyrrhenian Sea and North Africa). Temperature data was thus available for six of the study sites (Fig. 1) and only these study sites were included in the environmental association analyses. Of the temperature parameters collected from six of the sites (Table 3), no environmental associations were found with any loci.
Population genetic analysis using genome-wide SNPs identified from nine S. alveolata reef sites spanning the latitudinal and longitudinal range of the species revealed low gene flow between sites, equating to low effective dispersal, with 25 out of 36 pairwise FST values showing significance and an average pairwise FST of 0.28 ± 0.10 across the system. The DAPC supported the differentiation of the Tyrrhenian Sea site from the other sites, as also identified by high FST values between this site and all other sites (FST = 0.32 ± 0.12), and the STRUCTURE analysis revealed the North Atlantic site to be a separate population from all the others (FST = 0.42 ± 0.10). Genetic divergence in terms of FST was not significantly correlated with either straight-line distance or shortest ocean distance. Presence of low gene flow between reef sites is in line with a growing number of studies that have identified population structuring in marine species which, as with S. alveolata in our study, have a dispersing larval stage and no hard geographical barriers to dispersal [17, 18, 48,49,50]. For instance, Benestan et al. (2015)  identified strong evidence of weak neutral genetic structure (average pairwise FST = 0.0019) in the American lobster (H. americanus) along a sea surface temperature gradient using RAD genotyping despite a dispersing planktonic larval stage with a duration of 4–6 weeks. Similarly, the bat star Patiria miniate (Brandt, 1835) showed genetic population structure along a latitudinal gradient (average pairwise FST = 0.061) despite a longer larval dispersal period of 6–10 weeks , which is in line with the estimated larval duration of S. alveolata at 5–12 weeks, dependant on temperature and food availability [42, 51]. Indeed, a meta-analysis revealed that mean pelagic larval duration was not a good predictor of gene flow  and our findings support this interpretation.
Three sites were differentiated from all other sites: North Atlantic (average pairwise FST = 0.42 ± 0.10 and STRUCTURE), Tyrrhenian Sea (average pairwise FST = 0.32 ± 0.12 and DAPC) and Bay of Biscay (average pairwise FST = 0.25 ± 0.05). The observed differences in the results produced by the two cluster assignment methods can be explained by how they use data: STRUCTURE assigns cluster groupings in order to minimise Hardy–Weinberg and linkage disequilibria within clusters by looking at allele frequencies [53, 54]. In contrast, DAPC maximizes the variance in linear combinations of allele frequencies (principal components) between groups while minimizing variance within groups . Differences in genetic structure estimated by STRUCTURE and DAPC therefore likely result from how genetic variation is considered in either approach. The northern Irish Sea, southern Irish Sea (in terms of FST) and English Channel (in terms of FST and the DAPC) sites were also all differentiated from each other (Table 2; Fig. 3b). The ocean circulation modelling predicted high isolation (zero to low larval density; Additional file 1) among all release sites, even in the most optimistic (dispersive) scenario. Therefore, the FST results are supported, at least in part, by patterns of ocean circulation influencing passive larval dispersal along the Atlantic and Mediterranean coastlines of Europe and North Africa.
High population sub-structuring could be partially explained by two behavioural barriers to connectivity: first, spawning in polychaetes is triggered by the action of waves, especially during spring tides, which tends to push larvae toward the coast rather than offshore  and second, S. alveolata larvae act to reduce dispersal out of the bay where they were spawned by moving higher in the water column with an incoming- and lower with an outgoing-tide . This ‘tidal stream transport’, has been suggested as a possible mechanism for facilitating/restricting advective transport in a number of taxa (e.g. [56, 57]). Yet despite this, relatively high levels of genetic diversity were observed across the sites (He = 0.52 ± 0.04). Observed heterozygosity was very low at all sites and could be due to inbreeding , as a result of the site isolation predicted by the ocean circulation modelling and due to the reproductive strategy of gregarious colonial species, as juveniles settle together in patches and spawning of one individual triggers the spawning of those immediately adjacent . Heterozygote deficiency has also been recorded in marine invertebrates as a consequence of the Wahlund effect, due to the coexistence of genetically distinct cohorts within a sampling site . However, low heterozygosity estimations could also be a consequence of poor genome coverage and low sample size at some sites [59, 60]. Further research is needed to assess whether presence of null alleles is genome-wide, as evidence of inbreeding, or restricted to particular loci, as evidence of allelic dropout due to low coverage, in order to assess whether the species is at inbreeding risk.
Identifying the relative contribution of gene flow, genetic drift and natural selection to population structure is difficult in marine invertebrates due to their fluctuating population sizes . This is particularly difficult in S. alveolata, as in-depth local ecological knowledge, such as population size and breeding behaviour, is lacking in the majority of sites (but see [33, 61,62,63]). That said, the identification by STRUCTURE of only the North Atlantic site as a separate, but not completely isolated, population, with all other sites forming a single population, supports the idea that low levels of gene flow are still maintaining genetic diversity within the system as a whole. It could be that asymmetric dispersal, a common feature of seascapes due to unidirectional currents  has led to some sharing of genotypes between the North Atlantic and other sites, as observed in the STRUCTURE plot (Fig. 4b). Such wide-ranging genetic connectivity has also been observed in other polychaete species . This is potentially a positive sign for survival of the species in a changing climate, as the maintenance of genetic diversity is key to facilitating rapid evolution when environmental conditions change . However, S. alveolata reduce their larval duration with increasing temperatures  and our models show that shorter larval duration leads to reduced larval dispersal. Therefore, lower connectivity, and thus gene flow, is predicted for S. alveolata in a changing climate.
Despite some evidence of low levels of gene flow between the majority of the study sites within the system, both FST values and STRUCTURE analysis revealed a level of isolation of the North Atlantic site and identified this reef as a separate population to all others. There are two possible processes that could, separately or in synergy, be causing the observed pattern of population structure: historical isolation and contemporary patterns of gene flow. During the Last Glacial Maximum, geomorphology and fossil evidence suggests that southwest Ireland was partially unglaciated  and genetic data supports the presence of a glacial refugium in this area for both terrestrial [68, 69] and marine species [70,71,72]. In particular, Jolly et al. (2006)  found that two coastal polychaete worms (Lagis koreni Malmgren, 1866, formerly Pectinaria koreni, and Owenia fusiformis Delle Chiaje, 1844) showed a private Irish Sea haplotype linking two ancestral haplotypes and they suggest this could have evolved in a small ice-free area along the southwest coast of Ireland. Therefore, one potential explanation is that this population was isolated in a different glacial refugium to the rest of the sites, with subsequent admixture, but further study is needed to test this hypothesis.
The second hypothesis, that contemporary gene flow is very low between the North Atlantic and other S. alveolata sites, is supported by predictions of larval dispersal as seen in the ocean circulation modelling. In even the most optimistic scenario for larval dispersal, the North Atlantic site was not predicted to have interchange of individuals with any of the other sampled sites and had by far the highest predicted proportion of larval retention (79% of released individuals were retained compared to an average of 20 ± 13% for all other sites). This can be explained by the hydrodynamic modelling in terms of current patterns within Galway Bay limiting dispersal. The presence of the North Atlantic Current, which moves eastwardly towards Ireland and then continues northwards, is also an isolating factor for the North Atlantic reef. Any larvae that do move beyond their spawning reef at the North Atlantic site are drawn northwards, beyond the current northern range limit for the species (the Solway Firth, Scotland ). This oceanic barrier is also a likely cause of the observed genetic divergence between the northern/southern Irish Sea sites and the English Channel, and the reason that the dispersing larvae from more southerly sites do not reach the North Atlantic site. Although larval dispersal of many other marine species has been found to be mediated by the North Atlantic Current (e.g. [73, 74]), isolation due to the North Atlantic Current has not previously been reported, likely because it depends upon species-specific occurrence range and life history traits. Therefore, reduced contemporary gene flow is a likely cause of the observed population isolation of the North Atlantic S. alveolata reef. These findings highlight the North Atlantic reef site as at risk if population size reductions occur as recruitment and genetic augmentation is unlikely from elsewhere within the species’ range. Therefore, conservation management is needed to ensure that population size does not decrease at this site, and potentially throughout Ireland.
Despite the importance of ocean circulation to larval dispersal in general [18, 19, 75] and to S. alveolata in particular, as identified by our data in regards to the isolated North Atlantic site, our ocean circulation models of larval dispersal did not show a significant correlation with observed genetic distance between sites. This is surprising as both FST values and larval dispersal models suggest low connectivity between sites. This lack of correlation between our larval dispersal modelling and genetic distance values is in contrast to studies on a number of species with a dispersing larval stage including the Mediterranean shore crab (Carcinus aestuarii Nardo, 1847) , the bat star (P. miniata)  and the American lobster (H. americanus) ; these studies found genetic structure was directly related to ocean currents or to estimates of potential larval connectivity obtained with coupled physical-biological models. However, these studies were conducted over a smaller geographical area than our study. As per our study, Jorde et al. (2015)  found that when looking at large-scale differentiation patterns in the north Atlantic, geographic distance and larval drift alone explained only a minor portion (2.5–4.7%) of genetic isolation in the northern shrimp (Pandalus borealis Krøyer, 1838). Galindo et al. (2010)  modified a biophysical model for Monterey Bay in California to simulate dispersal of the acorn barnacle (Balanus glandula Darwin, 1854) but it also did not match an observed genetic cline in the species. They discovered that their model fit was improved by including natural selection, larval retention, and input values from an additional source population .
There are several factors likely to reduce the similarity between modelled larval dispersal and observed gene flow. Gene flow is representative of multigenerational mixing, while the model used in this study is representative of the dispersal patterns of only one generation. It is possible to simulate multigenerational gene flow based on dispersal models [49, 77, 78] but this would be confounded by the lack of known sites that could be used as “stepping stone” sites, which may facilitate dispersive spread over multiple generations. Such sites are likely to exist but are currently undocumented. These undiscovered sites represent missing source populations within our model and may well be a cause of the observed mismatch between predicted larval dispersal and observed gene flow. The distribution and occurrence of S. alveolata reefs are not well documented, particularly in the southern range of the species (Firth, unpublished), and the prediction of our model that genetic differentiation between our sample sites would be high could be due to the absence of intermediate reef sites within the model that would create higher levels of admixture within the system as a whole [9, 79]. This is further complicated by the fact that reefs can vary temporally in their presence at a site . Therefore, increased knowledge of the location of S. alveolata reefs and ongoing monitoring of reef sites is required to generate a clearer picture of the role of connectivity in meta-population maintenance in this species and to further inform the ocean circulation model. Our study highlights the importance of comparing bio-physical models with observed population structure of a species in order to create accurate dispersal predictions.
Loci under selection
Twenty-seven SNPs were identified as outlier loci, all of which were putatively under balancing selection. Identification of a higher number of loci potentially under balancing rather than divergent selection is common in genome scan studies, including those on marine invertebrates [12, 50, 80,81,82]. Our finding that 5.6% of SNPs were potentially under balancing selection is reflective of a growing body of evidence that suggests that balancing selection, which acts to preserve polymorphism, is more important in the genome than previously considered [12, 50, 83]. In their study of the sea anemone, Nematostella vectensis Stephenson, 1935, covering a broad geographical range, Reitzel et al. (2013)  identified 37 polymorphic sites inferred to be under balancing selection, but none under divergent selection, as with our study. However, further evidence in the form of elevated polymorphism, reduced differentiation, and shifts towards intermediate allele frequencies are needed to confirm that these loci are indeed under balancing selection .
The lack of identification of loci under divergent selection using either outlier analyses or environmental association tests is surprising given that local adaptation to temperature is known to be present in this system . We previously reported that S. alveolata individuals from sites along the latitudinal range of this species (all of which are included in this study) showed different responses to thermal regime changes in terms of their membrane lipid composition dependant on their site of origin . This is likely a reflection of poor genome coverage leading to a low number of genotyped loci (as per ) and whole genome sequencing would be beneficial to search for adaptive loci and their functions in this system .
Modelling of larval dispersal predicted low gene flow between sites for S. alveolata across their range, which was in part supported by the genomic data. In particular, population genetic analyses identified the North Atlantic reef is an isolated population, due to the effect on dispersal of the North Atlantic Current. These findings have implications for the management of S. alveolata engineered habitats in a changing climate and have highlighted the North Atlantic reef as at increased risk and in need of direct conservation action in order to preserve not only this species, but above all the crucial functional ecological roles (biodiversity hotspots, coastal protection, carbonate traps, etc.) that are associated with these bioconstructions. Our study highlights the utility of using seascape genomics to identify populations of conservation concern.
Sampling and molecular methods
In total, 180 mature individual worms were haphazardly collected from nine sites (20 per site) across the global geographic range between August 2013 and March 2014, from North Africa at the species’ southern boundary, along the Atlantic Coast, to the Northern Irish Sea in the north, as well as from two sites on the Mediterranean coast (Fig. 1). Individuals were placed immediately in RNA-later (Qiagen Inc., Crawley), maintained at room temperature for between 48 h and 1 week and then stored at − 20 °C until DNA extraction.
Genomic DNA was extracted using chaotropic “Chaos” buffer (4 M guanidine thiocyanate; 0.5% N-lauroyl sarcosine; 25 μM tris(hydroxymethyl)aminomethane (Tris), pH 8.0; 0.1 M 2-mercaptoethanol, 0.5% N-lauroyl sarcosine) (modified from Fukami et al. 2004 ). Worms were homogenised within the buffer at room temperature for at least 48 h before adding an equal volume of phenol extraction buffer (PEB) and double the volume of phenol:chloroform:isoamyl alcohol (25:24:1). From here the protocol followed a phenol-chloroform extraction with isopropanol precipitation . An additional ammonium acetate (3 M) precipitation step was added to further purify DNA. DNA was resuspended in sterile distilled H2O, quality assessed using a NanoDrop 8000 spectrophotometer (Thermo Scientific, Wilmington, DE), visualised on 1% agarose gels, and quantified using the Quant-IT PicoGreen dsDNA Assay Kit (Life Technologies, Carlsbad, California). Single-end Restriction-site Associated DNA (RAD) libraries were prepared by digesting DNA from each individual using the EcoRI and MseI restriction enzymes. Samples were individually barcoded and sequenced (Ecogenics, Zurich-Schlieren, Switzerland) using an Illumina NextSeq 500 (Illumina, San Diego, California).
SNP detection and summary statistics
All analyses were carried out on the Analyses and Bioinformatics for Marine Science (ABiMS) Galaxy Platform at the Station Biologique de Roscoff (galaxy.sb-roscoff.fr ;). Samples with mean Phred scores of < 25 or read coverage of < 150,000 were removed from the dataset and sequences were trimmed to 60 bp length using Trimmomatic version 0.36.4 . STACKS version 1.46.0  was used to de novo map the sequences and call SNPs by sequentially running “ustacks”, “cstacks” and “sstacks” using the parameters suggested to be optimal for population genetic inference by  namely, allowing a minimum depth of coverage of three (ustacks: m; stack coverage), a maximum of two mismatches between reads for a single individual (ustacks: M), and four mismatches between primary and secondary reads (ustacks: N). These settings have been found to increase the number of loci and reduce the SNP and allele error rates within a dataset . Furthermore, the low value of M used here reduces the risk of paralogous loci being merged into the same SNP . The number of mismatches allowed between loci was set at three (cstacks: n ;).The removal or break up of highly repetitive RAD tags was permitted within the program (ustacks: deleverage).
The programme populations from the STACKS suite, was used to process the SNP data . RADtags with at least 10x depth of read coverage were retained to ensure accuracy of heterozygous SNP calls (populations: p ;). Lumped paralogs can be a concern in high coverage data sets , but 10x coverage is considered low  and thus reduces this to a low risk factor. A locus had to be present in 70% of the populations to be included in the analysis, which allows for mutations in restriction sites that may cause loci to dropout in certain lineages [50, 94] and the minor allele frequency within populations (populations: min_maf) was set at > 0.01 . In order to maximise SNP discovery, the percentage of individuals required within a population to process a locus was set at 50% (populations: r). Only a single SNP per RADtag was retained for subsequent analyses to avoid problems of non-independence between markers . The STRUCTURE file output by populations was transformed in PGD Spider version 22.214.171.124  for use in all further analyses. Each population and locus was tested for deviation from Hardy-Weinberg equilibrium and linkage disequilibrium in ARLEQUIN version 126.96.36.199  and significance was assessed following Bonferroni correction for multiple tests. Loci that showed null alleles in multiple populations, or with significant deviation from Hardy–Weinberg equilibrium, were removed from further analyses to reduce the risk of inclusion of paralogous sequence variants .
Genomic summary statistics were calculated in ARLEQUIN version 188.8.131.52  as mean number of alleles per locus (N alleles), nucleotide diversity (θ and π), expected heterozygosity (He) and observed heterozygosity (Ho). Separate analyses of molecular variance (AMOVA) were carried out based on site of origin and genetic clusters identified in STRUCTURE (see below) using 16,000 permutations.
An outlier approach was used to identify loci that had a higher or lower FST than expected under a neutral model of selection (under positive or balancing selection, respectively). Outliers were estimated in two ways: firstly, using BAYESCAN version 2.1 , which uses a Bayesian approach to estimate the posterior probability that a locus is affected by selection, and was run using 20 pilot runs of 5000 iterations each, a total of 1,050,000 iterations (sample size of 100,000 and a thinning interval of 10), and a burn-in of 50,000. Only loci with a posterior probability (P) ≥ 0.95 with a prior odd of 10 were considered as outliers. BAYESCAN has consistently outperformed other outlier detection methods in terms of lack of false positives [4, 17, 98]. Outliers were also calculated in ARLEQUIN using 100,000 simulations, 500 demes, and a maximum expected heterozygosity of 0.5 . Outlier SNPs were identified using a P-value of ≤0.05. Only SNPs identified as outliers by both BAYESCAN and ARLEQUIN were conservatively selected as candidate loci under balancing or directional selection. Subsequent analyses of population structure and seascape genomics were conducted separately for neutral and outlier loci.
Pairwise FST were calculated in ARLEQUIN version 184.108.40.206  and significance in certainty of the estimator assessed following Bonferroni correction for multiple tests. Straight-line distances between sites were calculated using ArcGIS Version 10.1 (ESRI, 2011) and shortest ocean distances between sites were calculated in R version 3.1.2  using the package marmap , where distance was calculated excluding positive elevation . Isolation-by-distance (straight line and shortest ocean) was tested in ARLEQUIN (Version 220.127.116.11) using a Mantel test  with 10,000 permutations. To visualise these relationships, neighbor-joining (nj) tree estimation plots were constructed in R package ape 5.3 , following the methods outlined in , and were constructed using FST and shortest ocean distances for each pairwise (site) combination. Divergence between sites was also assessed using a Discriminant Analysis of Principal Components (DAPC) with the R package adegenet 2.1.2 , with 30 components and two discriminant functions retained in the analysis. Hierarchical structure was also considered by carrying out a second DAPC excluding the most divergent site, as identified in the first DAPC.
Genetic clusters (populations) were inferred using Bayesian analyses carried out in STRUCTURE version 2.3.4 . The analyses assumed admixture, correlated allele frequencies, and were run with 100,000 burn-in cycles and 100,000 Markov Chain Monte Carlo runs. The number of populations (K) was considered between 1 and 10, with 10 replicates per K. This was repeated with and without location priors, and using a hierarchical approach to resolve fine-scale genetic structure by removing highly differentiated populations . The most likely value of K was inferred by estimating ΔK using the Evanno method , as implemented in Structure Harvester version 0.6.94 .
Hydrodynamic models are increasingly used as predictors of larval dispersal (gene flow) in marine systems (e.g. [20, 107, 108] and many others). Probabilistic dispersal simulations of potential larval connectivity were run using the Connectivity Modeling System (CMS, ) paired with 3-hrly Hybrid Coordinate Ocean Model (HYCOM) and Navy coupled ocean data assimilation (NCODA) global 1/12° daily outputs from 2004, 2010, and 2012 . These years were selected to represent neutral, negative and positive phases of the North Atlantic Oscillation (NAO) respectively, and are thus likely to reflect average dispersal patterns over the last decade [109, 111]. Vertical velocities were calculated using the Eulerian continuity equation, and a random diffusive kick was added at each hourly dispersal time-step to better capture the potential effects of sub-gridscale turbulent diffusion (horizontal = 15 m2s− 1, vertical = 0.05 m2s− 1; after Kough et al. 2016 ). HYCOM is a global model appropriate for large landscape continental-scale dispersal simulations, but does not include tides. Larvae were considered as passive particles even though diel vertical migrations may occur . Simulations can therefore be over-dispersive, as the inclusion of tides and behaviour is liable to be retentive .
In the simulations, 100 propagules, each representing an S. alveolata larva, were released in a Lagrangian framework (so that each theoretical larvae was tracked individually), from each of nine locations (Fig. 1), and their dispersal under two different release scenarios predicted: (1) a conservative scenario with propagule releases every 6-h over 1 week in April for each year, centred on the spring tide, and (2) an optimistic scenario which modelled a daily release on every day of the year. The optimistic scenario was designed to simulate all the potential pathways of dispersal if reproduction occurs throughout the year, which has been suggested in S. alveolata , thus capturing the full variability of potential larval fates. Larval competency to settle was assumed after 4 weeks, and a maximum PLD of 12 weeks applied based on field and laboratory estimates [42, 51]. Larvae were permitted to settle when located within 5 km of a known reef site (Additional file 1). Additional (non-release) sites were also included as reported in The Global Biodiversity Information Facility (GBIF.org) in which larvae could also settle if encountered. Larval exchange counts (i.e. the source-sink contribution) were generated after 5, 8, 10, and 12 weeks PLD. Model outputs were processed in Matlab (Mathworks, Version R2016a). Pairwise matrices of larval exchange counts, standardised using proportions, were generated from larval fates averaged across the three focus years to compare against genetic pairwise FST matrices. Daily positions of larvae throughout the simulation were used to generate cumulative 2-D maps of larval density to highlight the predicted pathways of dispersal and areas of potential settlement.
The proportion of settlers reaching each study site from each release site under the four larval duration scenarios (5, 8, 10 and 12-weeks) in the ocean circulation modelling simulations were inverted and used to form a matrix of pairwise dispersal resistances (i.e. the proportion of larvae that did not reach each site). Mantel tests  were carried out in ARLEQUIN (Version 18.104.22.168)  using 16,000 permutations to assess the correlation between pairwise FST (using putatively neutral loci) and dispersal. A Bonferroni correction was applied to assess the significance of the Mantel tests .
Loci under selection
Locally collected temperature data provides insight into the conditions that S. alveolata experience at a spatial and temporal scale that cannot be achieved using satellite data  and can be used for comparison with genomic data, to test genotype-environment associations . Temperature data were obtained from Seabra et al. (2015)  who used biomimetic temperature loggers placed in the intertidal zone of exposed shores [115, 118]. At five sites, temperature data were collected with a resolution of 0.5 °C at mid-intertidal level every hour between July 2010 and July 2014. High- and low-tide temperatures were also identified to the nearest hour. At the remaining four sites (i.e. English Channel, Balearic Sea, Tyrrhenian Sea and North Africa), an iButton (Maxim, Munich) was placed in the mid-littoral zone, collecting temperature data with a resolution of 0.5 °C every 3 h between April 2014 and April 2015. Temperature readings closest to high- and low-tide were identified.
S. alveolata experience tidally-driven fluctuations in water and air temperatures. Sixteen metrics were calculated to describe variation in temperature at each site as follows: mean temperature, maximum temperature, minimum temperature, temperature range (maximum-minimum), standard deviation of the mean, 95th – 5th percentile (to exclude outlier temperatures), average daily temperature range (maximum-minimum per day), and average daily standard distribution (standard deviation of the mean), per site. High-tide (water) and low-tide (air) temperature values were also used separately to calculate mean temperature, maximum temperature, minimum temperature, and temperature range (maximum-minimum) for water and air temperatures, respectively (Table 3). Larvae of S. alveolata can be found in the water column all year round , suggesting that they are reproductively active throughout the year, thus temperature data from all months were included in the analyses.
We utilised a gene-environment association software, BAYENV2 , to test for covariance between SNP allele frequencies and the 16 thermal variables. BAYENV2 controls for neutral population structure by first estimating a covariance matrix among populations for all loci and then accounting for that covariance in the gene-environment association test [120, 121]. The programme was run with 100,000 iterations for both neutral parameterisation and association testing . Outlier loci were classed as those with a Bayes Factor > 3 [8, 123, 124].
Availability of data and materials
DNA sequences are available on NCBI SRA (http://www.ncbi.nlm.nih.gov/bioproject/645800; BioSample accessions: SAMN15516880–978) . Sample site locations, temperature data and ocean circulation model outputs are available on FIGSHARE (https://figshare.com/articles/dataset/Seascape_genomics_reveals_population_isolation_in_the_reef-building_honeycomb_worm_Sabellaria_alveolata_L_/12645896) . Galaxy workflow is available on github: https://github.com/annamuir/sabellaria.git
Planktonic larval duration
Restriction-site associated DNA sequencing
Reusch TBH, Wood TE. Molecular ecology of global change. Mol Ecol. 2007;16:3973–92.
Holderegger R, Kamm U, Gugerli F. Adaptive vs. neutral genetic diversity: implications for landscape genetics. Landsc Ecol. 2006;21:797–807.
Dionne M, Caron F, Dodson JJ, Bernatchez L. Landscape genetics and hierarchical genetic structure in Atlantic salmon: the interaction of gene flow and local adaptation. Mol Ecol. 2008;17:2382–96.
Chávez-Galarza J, Henriques D, Johnston JS, Azevedo JC, Patton JC, Muñoz I, De La Rúa P, Pinto MA. Signatures of selection in the Iberian honey bee (Apis mellifera iberiensis) revealed by a genome scan analysis of single nucleotide polymorphisms. Mol Ecol. 2013;22:5890–907.
Matala AP, Ackerman MW, Campbell MR, Narum SR. Relative contributions of neutral and non-neutral genetic differentiation to inform conservation of steelhead trout across highly variable landscapes. Evol Appl. 2014;7:682–701.
Conover DO, Schultz ET. Phenotypic similarity and the evolutionary significance of countergradient variation. Trends Ecol Evol. 1995;10:248–52.
Merilä J, Hendry AP. Climate change, adaptation, and phenotypic plasticity: the problem and the evidence. Evol Appl. 2014;7:1–14.
Pujolar JM, Jacobsen MW, Als TD, Frydenberg J, Munch K, Jõnsson B, Jian JB, Cheng L, Maes GE, Bernatchez L, Hansen MM. Genome-wide single-generation signatures of local selection in the panmictic European eel. Mol Ecol. 2014;23:2514–28.
Gleason LU, Burton RS. Genomic evidence for ecological divergence against a background of population homogeneity in the marine snail Chlorostoma funebralis. Mol Ecol. 2016;25:3557–73.
Lemay MA, Russello MA. Genetic evidence for ecological divergence in kokanee salmon. Mol Ecol. 2015;24:798–811.
Hohenlohe PA, Bassham S, Etter PD, Stiffler N, Johnson EA, Cresko WA. Population genomics of parallel adaptation in threespine stickleback using sequenced RAD tags. PLoS Genet. 2010;6:e1000862.
Reitzel AM, Herrera S, Layden MJ, Martindale MQ, Shank TM. Going where traditional markers have not gone before: utility of and promise for RAD sequencing in marine invertebrate phylogeography and population genomics. Mol Ecol. 2013;22:2953–70.
Sandoval-Castillo J, Robinson NA, Hart AM, Strain LWS, Beheregaray LB. Seascape genomics reveals adaptive divergence in a connected and commercially important mollusc, the greenlip abalone (Haliotis laevigata), along a longitudinal environmental gradient. Mol Ecol. 2018;27:1603–20.
Wernberg T, Smale DA, Thomsen MS. A decade of climate change experiments on marine organisms: procedures, patterns and problems. Glob Chang Biol. 2012;18:1491–8.
Munday PL, Warner RR, Monro K, Pandolfi JM, Marshall DJ. Predicting evolutionary responses to climate change in the sea. Ecol Lett. 2013;16:1488–500.
Reusch TBH. Climate change in the oceans: evolutionary versus phenotypically plastic responses of marine animals and plants. Evol Appl. 2014;7:104–22.
Benestan L, Quinn BK, Maaroufi H, Laporte M, Clark FK, Greenwood SJ, Rochette R, Bernatchez L, 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:5073–92.
Coscia I, Robins PE, Porter JS, Malham SK, Ironside JE. Modelled larval dispersal and measured gene flow: seascape genetics of the common cockle Cerastoderma edule in the southern Irish Sea. Conserv Genet. 2013;14:451–66.
Schiavina M, Marino IAM, Zane L, Melià P. Matching oceanography and genetics at the basin scale. Seascape connectivity of the Mediterranean shore crab in the Adriatic Sea. Mol Ecol. 2014;23:5496–507.
Xuereb A, Benestan L, Normandeau É, Daigle RM, Curtis JMR, Bernatchez L, Fortin MJ. Asymmetric oceanographic processes mediate connectivity and population genetic structure, as revealed by RADseq, in a highly dispersive marine invertebrate (Parastichopus californicus). Mol Ecol. 2018;27:2347–64.
Helmuth B, Mieszkowska N, Moore P, Hawkins SJ. Living on the edge of two changing worlds: forecasting the responses of rocky intertidal ecosystems to climate change. Annu Rev Ecol Evol Syst. 2006;37:373–404.
Harley CDG. Climate change, keystone predation, and biodiversity loss. Science. 2011;334:1124–7.
Marshall DJ, Monro K, Bode M, Keough MJ, Swearer S. Phenotype-environment mismatches reduce connectivity in the sea. Ecol Lett. 2010;13:128–40.
Jones C, Lawton J, Shachak M. Organisms as Ecosystem Engineers. In: Samson F, Knopf F, editors. Ecosystem Management. New York: Springer; 1996. p. 130–47.
Chaverra A, Wieters E, Foggo A, Knights AM. Removal of intertidal grazers by human harvesting leads to alteration of species interactions, community structure and resilience to climate change. Mar Environ Res. 2019;146:57–65.
Gruet Y. Spatio-temporal changes of sabellarian reefs built by the sedentary polychaete Sabellaria alveolata. Mar Ecol. 1986;7:303–19.
Wilson DP. The settlement behaviour of the larvae of Sabellaria alveolata (L.). J Mar Biol Assoc U K. 1968;48:387–435.
Wilson DP. Additional observations on larval growth and settlement of Sabellaria alveolata. J Mar Biol Assoc UK. 1970;50:1–31.
Buffet JP, Corre E, Duvernois-Berthet E, Fournier J, Lopez PJ. Adhesive gland transcriptomics uncovers a diversity of genes involved in glue formation in marine tube-building polychaetes. Acta Biomater. 2018;72:316–28.
Fournier J, Etienne S, Le Cam J-B. Inter- and intraspecific variability in the chemical composition of the mineral phase of cements from several tube-building polychaetes. Geobios. 2010;43:191–200.
Dubois SF, Commito JA, Olivier F, Retière C. Effects of epibionts on Sabellaria alveolata (L.) biogenic reefs and their associated fauna in the bay of Mont saint-Michel. Estuar Coast Shelf Sci. 2006;68:635–46.
Noernberg MA, Fournier J, Dubois SF, Populus J. Using airborne laser altimetry to estimate Sabellaria alveolata (Polychaeta: Sabellariidae) reefs volume in tidal flat environments Mauricio. Estuar Coast Shelf Sci. 2010;90:93–102.
Jones AG, Dubois SF, Desroy N, Fournier J. Interplay between abiotic factors and species assemblages mediated by the ecosystem engineer Sabellaria alveolata (Annelida: Polychaeta). Estuar Coast Shelf Sci. 2018;200:1–18.
Desroy N, Dubois SF, Fournier J, Ricquiers L, Le Mao P, Guerin L, Gerla D, Rougerie M, Legendre A. The conservation status of Sabellaria alveolata (L.) (Polychaeta: Sabellariidae) reefs in the bay of mont-Saint-Michel. Aquat Conserv Mar Freshwat Ecosyst. 2011;21:462–71.
Dubois SF, Barillé L, Cognie B. Feeding response of the polychaete Sabellaria alveolata (Sabellariidae) to changes in seston concentration. J Exp Mar Biol Ecol. 2009;376:94–101.
Porras R, Bataller JV, Murgui E, Torregrosa MT. Trophic structure and community composition of polychaetes inhabiting some Sabellaria alveolata (L.) reefs along the Valencia Gulf Coast, Western Mediterranean. Mar Ecol. 1996;17:583–602.
Dubois SF, Retière C, Olivier F. Biodiversity associated with Sabellaria alveolata (Polychaeta: Sabellariidae) reefs: effects of human disturbances. J Mar Biol Assoc U K. 2002;82:817–26.
De Smet B, D’Hondt AS, Verhelst P, Fournier J, Godet L, Desroy N, Rabaut M, Vincx M, Vanaverbeke J. Biogenic reefs affect multiple components of intertidal soft-bottom benthic assemblages: the Lanice conchilega case study. Estuar Coast Shelf Sci. 2015;152:44–55.
Farcy S. Analyse de la structure génétique des récifs de Sabellaria alveolata (L.) dans la baie du Mont-Saint-Michel à l’aide de marqueurs microsatellites (masters thesis): Université Pierre et Marie Curie – Paris VI; 2003.
Rigal F. Barrières biogéographiques et processus historiques chez les invertébrés marins : définition des unités taxonomiques et populationnelles chez Sabellaria alveolata (masters thesis): Université de Paris XI; 2005.
Pawlik JR. Larval settlement and metamorphosis of two gregarious Sabellariid polychaetes: Sabellaria alveolata compared with Phragmatopoma californica. J Mar Biol Assoc U K. 1988;68:101–24.
Dubois SF, Comtet T, Retière C, Thiébaut E. Distribution and retention of Sabellaria alveolata larvae (Polychaeta: Sabellariidae) in the bay of Mont-saint-Michel. France. Mar Ecol Prog Ser. 2007;346:243–54.
Muir AP, Nunes FLD, Dubois SF, Pernet F. Lipid remodelling in the reef-building honeycomb worm, Sabellaria alveolata, reflects acclimation and local adaptation to temperature. Sci Rep. 2016;6:1–10.
Team RC. R: A Language and Environment for Statistical Computing. In: R Foundation for Statistical Computing Vienna, Austria; 2014. https://www.R-project.org.
Muir AP, Thomas R, Biek R, Mable BK. Using genetic variation to infer associations with climate in the common frog, Rana temporaria. Mol Ecol. 2013;22:3737–51.
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:977–93.
Robins PE, Neill SP, Giménez L, Stuart R, Jenkins SR, Malham SK. Physical and biological controls on larval dispersal and connectivity in a highly energetic shelf sea. Limnol Oceanogr. 2013;58:505–24.
Galindo HM, Pfeiffer-Herbert AS, McManus MA, Chao Y, Chai F, Palumbi SR. Seascape genetics along a steep cline: using genetic patterns to test predictions of marine larval dispersal. Mol Ecol. 2010;19:3692–707.
Sunday JM, Popovic I, Palen WJ, Foreman MGG, Hart MW. Ocean circulation model predicts high genetic structure observed in a long-lived pelagic developer. Mol Ecol. 2014;23:5036–47.
Benestan L, Gosselin T, Perrier C, Sainte-Marie B, Rochette R, Bernatchez L. RAD genotyping reveals fine-scale genetic structuring and provides powerful population assignment in a widely distributed marine species, the American lobster (Homarus americanus). Mol Ecol. 2015;24:3299–315.
Cazaux C. Developement larvaire de Sabellaria alveolata. Oceanography Monaco. 1964;62:1–15.
Weersing K, Toonen RJ. Population genetics, larval dispersal, and connectivity in marine systems. Mar Ecol Prog Ser. 2009;393:1–12.
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:94.
Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genet Mol Res. 2000;155:945–59.
McCarthy DA, Young CM, Emson RH. Influence of wave-induced disturbance on seasonal spawning patterns in the sabellariid polychaete Phragmatopoma lapidosa. Mar Ecol Prog Ser. 2003;256:123–33.
Tankersley RA, Forward RB Jr. Endogenous swimming rhythms in estuarine crab megalopae: implications for flood-tide transport. Mar Biol. 1994;118:415–23.
Knights AM, Crowe TP, Burnell G. Mechanisms of larval transport: vertical distribution of bivalve larvae varies with tidal conditions. Mar Ecol Prog Ser. 2006;326:167–74.
Buckley J, Holub EB, Koch MA, Vergeer P, Mable BK. Restriction associated DNA-genotyping at multiple spatial scales in Arabidopsis lyrata reveals signatures of pathogen-mediated selection. BMC Genomics. 2018;19:1–21.
Benestan L, Ferchaud AL, Hohenlohe PA, Garner BA, Naylor GJP, Baums IB, Schwartz MK, Kelley JL, Luikart G. Conservation genomics of natural and managed populations: building a conceptual and practical framework. Mol Ecol. 2016;25:2967–77.
Nielsen R, Paul JS, Albrechtsen A, Song YS. Genotype and SNP calling from next-generation sequencing data. Nat Rev Genet. 2011;12:443–51.
Wilson DP. Sabellaria colonies at duckpool, North Cornwall, 1961–1970. J Mar Biol Assoc U K. 1971;51:509–80.
Firth LB, Mieszkowska N, Grant LM, Bush LE, Davies AJ, Frost MT, Moschella PS, Burrows MT, Cunningham PN, Dye SR, Hawkins SJ. Historical comparisons reveal multiple drivers of decadal change of an ecosystem engineer at the range edge. Ecol Evol. 2015;5:3210–22.
Curd A, Pernet F, Corporeau C, Delisle L, Firth LB, Nunes FLD, Dubois SF. Connecting organic to mineral: how the physiological state of an ecosystem-engineer is linked to its habitat structure. Ecol Indic. 2019;98:49–60.
Buonomo R, Assis J, Fernandes F, Engelen AH, Airoldi L, Serrão EA. Habitat continuity and stepping-stone oceanographic distances explain population genetic connectivity of the brown alga Cystoseira amentacea. Mol Ecol. 2017;26:766–80.
Nunes FLD, Van Wormhoudt A, Faroni-Perez L, Fournier J. Phylogeography of the reef-building polychaetes of the genus Phragmatopoma in the western Atlantic region. J Biogeogr. 2017;44:1612–25.
Mergeay J, Santamaría L. Evolution and biodiversity: the evolutionary basis of biodiversity and its potential for adaptation to global change. Evol Appl. 2012;5:103–6.
Maggs CA, Castilho R, Foltz D, Henzler C, Jolly MT, Kelly J, Olsen J, Perez K, Stam W, Vainola R, Viard F, Wares J. Evaluating signatures of glacial refugia for North Atlantic marine organisms. Ecology. 2008;89:S108–22.
Rowe G, Harris DJ, Beebee TJC. Lusitania revisited: a phylogeographic analysis of the natterjack toad Bufo calamita across its entire biogeographical range. Mol Phylogenet Evol. 2006;39:335–46.
Teacher AGF, Garner TWJ, Nichols RA. European phylogeography of the common frog (Rana temporaria): routes of postglacial colonization into the British Isles, and evidence for an Irish glacial refugium. Heredity. 2009;102:490–6.
Jolly MT, Viard F, Gentil F, Thiébaut E, Jollivet D. Comparative phylogeography of two coastal polychaete tubeworms in the Northeast Atlantic supports shared history and vicariant events. Mol Ecol. 2006;15:1841–55.
Hoarau G, Coyer JA, Veldsink JH, Stam WT, Olsen JL. Glacial refugia and recolonization pathways in the brown seaweed Fucus serratus. Mol Ecol. 2007;16:3606–16.
Provan J, Bennett KD. Phylogeographic insights into cryptic glacial refugia. Trends Ecol Evol. 2008;23:564–71.
Scheltema R. Larval dispersal as a means of genetic exchange between geographically separated populations of shallow-water benthic marine gastropods. Biol Bull. 1971;140:284–322.
Kettle AJ, Haines K. How does the European eel (Anguilla anguilla) retain its population structure during its larval migration across the North Atlantic Ocean? Can J Fish Aquat Sci. 2006;63:90–106.
Cowen RK, Sponaugle S. Larval dispersal and marine population connectivity. Annu Rev Mar Sci. 2009;1:443–66.
Jorde PE, Søvik G, Westgaard JI, Albretsen J, André C, Hvingel C, Johansen T, Sandvik AD, Kingsley M, Jørstad KE. Genetically distinct populations of northern shrimp, Pandalus borealis, in the North Atlantic: adaptation to different temperatures as an isolation factor. Mol Ecol. 2015;24:1742–57.
Kool JT, Paris CB, Andréfouët S, Cowen RK. Complex migration and the development of genetic structure in subdivided populations: an example from Caribbean coral reef ecosystems. Ecography. 2010;33:597–606.
Foster NL, Paris CB, Kool JT, Baums IB, Stevens JR, Sanchez JA, Bastidas C, Agudelo C, Bush P, Day O, Ferrari R, Gonzalez P, Gore S, Guppy R, McCartney MA, McCoy C, Mendes J, Srinivasan A, Steiner S, Vermeij MJA, Weil E, Mumby PJ. Connectivity of Caribbean coral populations: complementary insights from empirical and modelled gene flow. Mol Ecol. 2012;21:1143–57.
Kimura M, Weiss GH. The stepping stone model of population structure and the decrease of genetic correlation with distance. Genetics. 1964;49:561–76.
Defaveri J, Jonsson PR, Merilä J. Heterogeneous genomic differentiation in marine threespine sticklebacks: adaptation along an environmental gradient. Evolution. 2013;67:2530–46.
Vera M, Díez-Del-Molino D, García-Marín JL. Genomic survey provides insights into the evolutionary changes that occurred during European expansion of the invasive mosquitofish (Gambusia holbrooki). Mol Ecol. 2016;25:1089–105.
Moore JS, Bourret V, Dionne M, Bradbury I, O’Reilly P, Kent M, Chaput G, Bernatchez L. Conservation genomics of anadromous Atlantic salmon across its North American range: outlier loci identify the same patterns of population structure as neutral loci. Mol Ecol. 2014;23:5680–97.
Charlesworth D. Balancing selection and its effects on sequences in nearby genome regions. PLoS Genet. 2006;2:379–84.
Mable BK. Conservation of adaptive potential and functional diversity: integrating old and new approaches. Conserv Genet. 2019;20:89–100.
Fukami HH, Budd AFAF, Levitan DRDR, Jara JJ, Kersanach RR, Knowlton NN. Geographic differences in species boundaries among members of the Montastraea annularis complex based on molecular and morphological markers. Evolution. 2004;58:324–37.
Nunes FLD, Fukami HH, Vollmer SV, Norris RD, Knowlton N. Re-evaluation of the systematics of the endemic corals of Brazil by molecular data. Coral Reefs. 2008;27:423–32.
Le Bras Y, Roult A, Monjeaud C, Bahin M, Quénez O, Hériveau C, Bretaudeau A, Sallou O, Collin O. Towards a life sciences virtual research environment. An e-Science initiative in Western France. JOBIM. 2013:1–10.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Catchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA. Stacks: an analysis tool set for population genomics. Mol Ecol. 2013;22:3124–40.
Mastretta-Yanes A, Arrigo N, Alvarez N, Jorgensen TH, Piñero D, Emerson BC. Restriction site-associated DNA sequencing, genotyping error estimation and de novo assembly optimization for population genetic inference. Mol Ecol Resour. 2015;15:28–41.
Miller AD, van Rooyen A, Rašić G, Ierodiaconou DA, Gorfine HK, Day R, Wong C, Hoffmann AA, Weeks AR. Contrasting patterns of population connectivity between regions in a commercially important mollusc Haliotis rubra: integrating population genetics, genomics and marine LiDAR data. Mol Ecol. 2016;25:3845–64.
Matz MV. Fantastic beasts and how to sequence them: ecological genomics for obscure model organisms. Trends Genet. 2017;34:121–32.
Paris JR, Stevens JR, Catchen JM. Lost in parameter space: a road map for stacks. Methods Ecol Evol. 2017;8:1360–73.
Jeffries DL, Copp GH, Handley LL, Håkan Olsén K, Sayer CD, Hänfling B. Comparing RADseq and microsatellites to infer complex phylogeographic patterns, an empirical perspective in the Crucian carp, Carassius carassius. L. Mol Ecol. 2016;25:2997–3018.
Lischer HEL, Excoffier L. PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics. 2012;28:298–9.
Excoffier L, Lischer HEL. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and windows. Mol Ecol Resour. 2010;10:564–7.
Campbell NR, LaPatra SE, Overturf K, Towner R, Narum SR. Association mapping of disease resistance traits in rainbow trout using restriction site associated DNA sequencing. G3. 2014;4:2473–81.
Narum SR, Hess JE. Comparison of FST outlier tests for SNP loci under selection. Mol Ecol Resour. 2011;11(SUPPL. 1):184–94.
Van Wyngaarden M, Snelgrove PVR, DiBacco C, Hamilton LC, Rodríguez‐Ezpeleta N, Jeffery NW, RRE S, Bradbury IR. Identifying patterns of dispersal, connectivity and selection in the sea scallop, Placopecten magellanicus, using RADseq-derived SNPs. Evol Appl. 2017;10:102–17.
Pante E, Simon-Bouhet B. marmap: A package for importing, plotting and analyzing bathymetric and topographic data in R. PLoS One. 2013;8:6–9.
Mantel N, Valand R. A technique of nonparametric multivariate analysis. Biometrics. 1970;26:547–58.
Paradis E, Schliep K. Ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics. 2019;35:526–8.
Saitou N, Nei M. The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol. 1987;4:406–25.
Jombart T. Adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24:1403–5.
Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14:2611–20.
Earl DA. vonHoldt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012;4:359–61.
North EW, Schlag Z, Hood RR, Li M, Zhong L, Gross T, Kennedy VS. Vertical swimming behavior influences the dispersal of simulated oyster larvae in a coupled particle-tracking and hydrodynamic model of Chesapeake Bay. Mar Ecol Prog Ser. 2008;359(Leis 2007):99–115.
Truelove NK, Box SJ, Aiken KA, Blythe-Mallett A, Boman EM, Booker CJ, Byfield TT, Cox CE, Davis MH, Delgado GA, Glazer BA, Griffiths SM, Kitson-Walters K, Kough AS, Pérez Enríquez R, Preziosi RF, Roy ME, Segura-García I, Webber MK, Stoner AW. Isolation by oceanic distance and spatial genetic structure in an overharvested international fishery. Divers Distrib. 2017;23:1292–300.
Fox AD, Henry LA, Corne DW, Roberts JM. Sensitivity of marine protected area network connectivity to atmospheric variability. R Soc Open Sci. 2016;3:160494.
Chassignet EP, Hurlburt HE, Smedstad OM, Halliwell GR, Hogan PJ, Wallcraft AJ, Baraille R, Bleck R. The HYCOM (HYbrid Coordinate Ocean model) data assimilative system. J Mar Syst. 2007;65:60–83.
Ross RE, Nimmo-Smith WAM, Howell KL. Increasing the depth of current understanding: sensitivity testing of deep-sea larval dispersal models for ecologists. PLoS One. 2016;11:1–25.
Kough A, Claro R, Lindeman K, Paris C. Decadal analysis of larval connectivity from Cuban snapper (Lutjanidae) spawning aggregations based on biophysical modeling. Mar Ecol Prog Ser. 2016;550:175–90.
Müller M, Haak H, Jungclaus JH, Sündermann J, Thomas M. The effect of ocean tides on a climate model simulation. Ocean Model. 2010;35:304–13.
Muir AP, Biek R, Thomas R, Mable BK. Local adaptation with high gene flow: temperature parameters drive adaptation to altitude in the common frog (Rana temporaria). Mol Ecol. 2014;23:561–74.
Seabra R, Wethey DS, Santos AM, Lima FP. Side matters: microhabitat influence on intertidal heat stress over a large geographical scale. J Exp Mar Biol Ecol. 2011;400:200–8.
Bernatchez S, Xuereb A, Laporte M, Benestan L, Steeves R, Laflamme M, Bernatchez L, Mallet MA. Seascape genomics of eastern oyster (Crassostrea virginica) along the Atlantic coast of Canada. Evol Appl. 2019;12:587–609.
Seabra R, Wethey DS, Santos AM, Lima FP. Understanding complex biogeographic responses to climate change. Sci Rep. 2015;5:3–8.
Lima FP, Wethey DS. Robolimpets: measuring intertidal body temperatures using biomimetic loggers. Limnol Oceanogr Methods. 2009;7:347–53.
Coop G, Witonsky D, Di Rienzo A, Pritchard JK. Using environmental correlations to identify loci underlying local adaptation. Genetics. 2010;185:1411–23.
Bernatchez S, Laporte M, Perrier C, Sirois P, Bernatchez L. Investigating genomic and phenotypic parallelism between piscivorous and planktivorous lake trout (Salvelinus namaycush) ecotypes by means of RADseq and morphometrics analyses. Mol Ecol. 2016;25:4773–92.
Lotterhos KE, Whitlock MC. The relative power of genome scans to detect local adaptation depends on sampling design and statistical method. Mol Ecol. 2015;24:1031–46.
Wenzel MA, Piertney SB. Fine-scale population epigenetic structure in relation to gastrointestinal parasite load in red grouse (Lagopus lagopus scotica). Mol Ecol. 2014;23:4256–73.
Schweizer RM, VonHoldt BM, Harrigan R, Knowles JC, Musiani M, Coltman D, Novembre J, Wayne RK. Genetic subdivision and candidate genes under selection in North American grey wolves. Mol Ecol. 2016;25:380–402.
De La Torre AR, Roberts DR, Aitken SN. Genome-wide admixture and ecological niche modelling reveal the maintenance of species boundaries despite long history of interspecific gene flow. Mol Ecol. 2014;23:2046–59.
Muir AP, Dubois SF, Ross RE, Firth LB, Knights AM, Lima FP, Seabra R, Corre E, Le Corguillé G, Nunes FLD. Data from: Seascape genomics reveals population isolation in the reef-building honeycomb worm, Sabellaria alveolata (L.). NHBS SRA. 2020; https://www.ncbi.nlm.nih.gov/bioproject/645800.
Muir AP, Dubois SF, Ross RE, Firth LB, Knights AM, Lima FP, Seabra R, Corre E, Le Corguillé G, Nunes FLD. Data from: Seascape genomics reveals population isolation in the reef-building honeycomb worm, Sabellaria alveolata (L.). Figshare. 2020; https://figshare.com/articles/dataset/Seascape_genomics_reveals_population_isolation_in_the_reef-building_honeycomb_worm_Sabellaria_alveolata_L_/12645896.
We thank Yannis Moreau, Hugo Covés, Victoria Paterson, Francesca Colombo, Aurélie Barraud, Paulo Chainho, Cédric Hennache, Jose Luis Acuna, Nina Schlapfes, and Barbara La Porta for their help with field sampling. Thanks to Nataliya Stashchuk for useful discussions about the ocean circulation modelling and to Barbara Mable for useful comments on the manuscript. Sampling in Morocco was carried out in collaboration with Hakima Zidane and Chlaida Malika of the National Fisheries Research Institute, Casablanca, Morocco.
Fieldwork, laboratory work, and genomic analyses: A. Muir and F. Nunes (FN) were supported by a grant from the Regional Council of Brittany, from the European Regional Development Fund (ERDF) and supported by the “Laboratoire d’Excellence” LabexMER (ANR-10-LABX-19) and co-funded by a grant from the French government under the program “Investissements d’Avenir”, to FN. Ocean circulation modelling: R. Ross and A. Knights (AK), were supported by an internal grant from the Faculty of Science and Environment, University of Plymouth, to AK. Temperature data collection: F. P. Lima and R. Seabra were supported by FEDER funds through COMPETE (POCI-01-0145-FEDER-031088) and through POR Norte (NORTE-01-0145-FEDER-031053), and by Portuguese national funds through FCT (PTDC/BIA-BMA/31088/2017, PTDC/BIA-BMA/31053/2017). The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Ethics approval and consent to participate
Permissions were not required to collect the samples in this study.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Muir, A.P., Dubois, S.F., Ross, R.E. et al. Seascape genomics reveals population isolation in the reef-building honeycomb worm, Sabellaria alveolata (L.). BMC Evol Biol 20, 100 (2020). https://doi.org/10.1186/s12862-020-01658-9
- Ocean circulation modelling
- Marine invertebrate
- Larval dispersal