Recent papers postulate that epifaunal organisms use artificial structures as stepping-stones to spread to areas that are too distant to reach in a single generation. With thousands of artificial structures present in the North Sea, we test the hypothesis that these structures are connected by water currents and act as an interconnected reef. Population genetic structure of the Blue mussel, Mytilus edulis was expected to follow a pattern predicted by particle tracking models (PTM). Correlation between population genetic differentiation, based on microsatellite markers, and particle exchange was tested. Specimens of M. edulis were found at each location, although the PTM indicated that locations >85 km offshore were isolated from coastal sub-populations. Fixation coefficient FST correlated with the number of arrivals in the PTM. However, the number of effective migrants per generation as inferred from coalescent simulations did not show a strong correlation with the arriving particles. Isolation by distance analysis showed no increase in isolation with increasing distance and we did not find clear structure among the populations. The marine stepping-stone effect is obviously important for the distribution of M. edulis in the North Sea and it may influence ecologically comparable species in a similar way. In the absence of artificial shallow hard substrates, M. edulis would be unlikely to survive in offshore North Sea waters. Although we found an indication that FST was lower between connected locations, isolation by distance analysis showed no increase in isolation with increasing distance. Finally, we did not find clear structure among the populations.

Mytilus edulis population genetics

Samples were collected at 27 locations (henceforth named sub-populations) between April 2014 and January 2016: from coastal sub-populations during low tide, by commercial and scientific divers from oil and gas platforms and wind farm foundations and during inspection, repair and maintenance work on buoys and offshore installations from structures lifted out of the water (Table 1, Figure 1). Sampling depth varied between 0 and 27 meters. Two likely outgroup sub-populations corresponding to peripheral populations of the species M. edulis (Denmark) and M. galloprovincialis (Portugal) were included; from a harbour in Lisbon, Portugal and from mussel longlines in the Limfjorden in Denmark. Geographic distances between sub-populations (excluding outgroups) varied from 17 km (Scheveningen and Q13-A) to 1,105 km (Blyth and Sylt) with an average of 348 km.

Between 50 and approximately 100 individuals were sampled randomly at every location and stored at -20°C or in 70% ethanol to be transported to the laboratory. Then, all samples were cleaned of marine growth and placed at -80°C for long-term storage.

Molecular methods

From each sample, between 24 and 67 specimens were selected randomly and genomic DNA was isolated from the adductor muscles by using the 96 well genomic DNA extraction kit according to the manufacturer’s protocol (FAVORGEN Biotech Corp). DNA concentration was measured on a Nanodrop and diluted to 10 ng/μl. Two multiplex sets of four markers were used (set1: Med367, Med379, Med722 and Med733; set2: Med737, Med740, Med747 and Me15/16). Markers were used individually in the PCR, pooled per set and analysed on an ABI3730 DNA Analyser. Seven markers were microsatellite loci as described in Lallias et al. (2009) and the Me15/16 locus targeted a part of the adhesive protein gene which partially discriminates between M. edulis, M. galloprovincialis and M. trossulus (Inoue, Waite, Matsuoka, Odo, & Harayama, 1995). Marker information, PCR conditions and multiplex conditions are indicated in Table 2. The GeneScan 500 LIZ marker was used as internal marker. Allele calling was performed in Genemapper v3.7 (Applied Biosystems 2004).

Connectivity calculated by particle tracking models

The transport of mussel larvae between sampled sub-populations in the North Sea was modelled using two Delft3D software modules: FLOW and PART (Deltares, 2016b, 2016a). The PART module was able to simulate mid-field water quality and particle tracking, based on a hydrodynamic forcing output from the other Delft3D module, FLOW.

Hydrodynamical model

The hydrodynamical model applied was Delft3D-FLOW, which solves the unsteady shallow-water equations in three dimensions. The model incorporates a large number of processes, such as wind shear, wave forces, tidal forces, density-driven flows, stratification, atmospheric pressure changes, air temperature and the exposure and inundation of intertidal flats. The large number of processes included in this module means that Delft3D-FLOW can be applied to a wide range of environments (e.g. river, estuarine, coastal and marine areas; Lesser et al. 2004). Flow equations were solved on a curvilinear grid consisting of 8,710 computational elements (Roelvink, Jeuken, van Holland, Aarninkhof, & Stam, 2001). The vertical resolution of the model was ten water layers using a sigma-coordinated approach (i.e., proportional to water depth; Stelling and van Kester 1994). Hydrodynamic transport was computed using detailed bathymetry and open boundary forcing based on tidal constituents. The model was forced using meteorological data from the High Resolution Limited Area Model (KNMI, 2015), which comprised two horizontal wind velocity components (at 10 m above mean sea level) and other atmospheric variables such as air pressure and temperature, archived every 6 h. The freshwater discharges from 18 rivers were included in the model; seven of these discharges varied temporally (daily averages) and 11 were constant (based on long-term averages). This model is described in detail in Erftemeijer et al. (2009).

Two grid lay-outs covering the southern North Sea (including the Wadden Sea) were used in this study: a moderately fine grid (ZUNOGROF) and a domain decomposition model grid (ZUNO-DD) with a much higher grid resolution in the Dutch coastal zone and the Wadden Sea. The subdomains of the ZUNO DD model are displayed in Figure 2 using different colours and are used to increase spatial resolution for hydrodynamic results in these areas. In coastal and inshore (estuarine, lagoonal) systems, the spatial variability in hydrodynamic forcing is much higher. To enable the simulation of these processes, the areas near the Dutch coast and in the Wadden Sea have been given a much higher horizontal grid density, so with smaller grid cells per surface unit. This provides for a much better simulation of the hydrodynamic processes in such areas. These spatial differences in resolution in the Dutch coastal waters and the Wadden Sea have been given a different colouring in Figure 2.

Particle transport model

Particle tracking models (PTM) are often used in environmental modelling (e.g., North et al. 2008; Broekhuizen et al. 2011; Postma et al. 2013). Here, the Delft3D-PART module was used to calculate larval transport across the southern North Sea. Delft3D-PART is a random walk particle tracking model, based on the principle that movement of dissolved substances in water can be described by a limited but potentially large number of discrete particles that are subject to advection due to currents and by horizontal and vertical dispersion. The movement of the particles in the model consists of two steps: advection, which is driven by the FLOW results, and dispersion, which is a stochastic random walk process. In addition, the horizontal and vertical movement of the particles can be adjusted to account for swimming preferences or changes in buoyancy. Particle tracking allows water quality processes to be described in a detailed spatial pattern, resolving sub-grid concentration distributions. Delft3D-PART is shown to be locally mass conservative (Postma et al., 2013).

Modelling set up

DELFT3D-FLOW calculated the hydrodynamic conditions for the southern North Sea for two consecutive years, 2004 and 2005. Meteorological conditions are important to determine the geographic destination of particles such as larvae, superimposed on the regular transport conditions in the southern North Sea. At the time of analysis, there were no hydrodynamic modelling results available for the years of sampling, 2014 to 2016. Instead, two model runs from 2004 and 2005 were used. Although there is large variation in climatic forcing between years, we chose these two years because of the overall comparable (southern) North Sea wind patterns during the months we used for PTM, i.e. March to June. The comparability especially regards wind directions and wind force. We could only roughly assess this comparability from meteorological data available to us (Royal Dutch Meteorological Institute); we did not do a detailed analysis. We assume that through this approach the differences between the hydrodynamic circulation patterns in the (southern) North Sea between these two periods is less than when random years would have been used. For the DELFT3D-PART model, each sampling location acted as an origin point of larvae and the density of the larvae throughout the consecutive months at the other locations (destinations) was calculated. Larval release was started at day 59 (1 March) of each year, from the top layer in the model. Both the timing and the period of particle release was based on the general reproductive timing and pelagic larval stage period of M. edulis. Mussels start spawning early spring, based on temperature of the water. Variation of spawning timing may result of the variation in temperature increase between years; in the southern North Sea this usually starts in March, and peaks in April/May or even later (De Vooys, 1999). Time from spawning of mussel larvae to settlement is around three months (Widdows, 1991). It was further assumed that the larvae are largely passive and that they are transported only by advection and general dispersion. They were also considered neutrally buoyant. Each origin released 1,489 particles per half hour, for 14 days (1,000,001 particles in total), after which no further release of particles occurred. The program calculated the horizontal transport and vertical dispersion to lower layers with time steps of 30 minutes for 70 days. Although mortality of M. edulis larvae is expected to be very high (Barker Jørgensen, 1981), this study did not include this aspect in the DELFT3D-PART model. The purpose of the study was to assess the maximum success of larvae settling at a specific destination, which was assumed to be proportional to the total number of larvae arriving at a specific destination throughout the 70 days. The total number of arriving particles from each source-destination combination was calculated and combined in a particle matrix.

Data analysis

Hybrid removal

A mix of alleles of Mytilus galloprovincialis and M. trossulus were observed in our data at the partially species diagnostic locus Me15/16, and therefore we estimated the hybrid index h for each individual using a maximum likelihood approach that estimates the proportion of alleles inherited from one of two hybridizing species (Buerkle, 2005; Buerkle & Lexer, 2008). This approach was applied as implemented in the R package Introgress (Gompert & Alex Buerkle, 2010). The samples from Lisbon (M. galloprovincialis) and Limfjorden (M. edulis) were used as parental reference samples and locus Med379 was not taken into account due to high levels of missing data for several samples. All individuals identified as non-pure M. edulis (h < 1.0) were omitted from subsequent analyses, as this aspect would have affected the non-neutral influences on connectivity inferences.

Isolation with migration analysis

Effective gene flow between sampling locations was estimated for all possible pairs, using the coalescent approach implemented in the software 'Isolation with Migration' (IMa2p) (Hey & Nielsen, 2004; Sethuraman & Hey, 2016). The IMa2p model consists of subpopulations that are simulated backwards in time towards merging sometime in the past. Parameters for time since population divergence t, population sizes q and bidirectional migration rates m are estimated by a Metropolis-coupled Markov Chain Monte Carlo process. In order to not overfit the data we only ran simulations for pairs of samples. Coalescent simulation runs consisted of ten Markov Chain Monte Carlo chains with geometric heating (setting parameters ha = 0.99, hb = 0.75) and ten million steps after an initial burn-in period of five million steps. Parameters estimated and parameter ranges examined (in coalescent units) were three population sizes q (ancestral and two populations ensuing from population subdivision, range examined 0-100), two migration rates after population subdivision (both directions, range examined 0-1) and time since population subdivision t (range examined 0-10). Convergence of estimated parameter distributions was ensured by: checking effective sample size (ESS) values, autocorrelation values and chain swapping, examining trend line plots for absence of trends, and by comparing parameter estimates generated from the genealogies produced during the first and second half of runs. A mutation rate for the seven microsatellite loci of 0.001564 per generation was assumed, following the methodology by Luttikhuizen et al. (2018). Migration rates among sampled locations were visualised by chord diagrams using the R package Circlize version 0.4.3 (Gu, Gu, Eils, Schlesner, & Brors, 2014).

The fit of the pairwise isolation-with-migration models to the data and the underlying processes was assessed by comparing the distribution of estimated migration rates with that of the pairwise F_ST values directly estimated from the microsatellite data. Pairwise F_ST values were converted to 2N_em values following the infinite island model approximation by Wright (Wright, 1943) in which F_ST = (1 + 4N_em)^-1. Correspondence between 2N_em values stemming from F_ST estimation and from isolation-with-migration simulations was compared using Pearson rank correlation. Correlation was performed with the average between the two (bidirectional) migration rates stemming from isolation-with-migration simulations.

EXPLANATION OF DATASET
The data are formatted in csv files contained in two zip files:
Coolen_et_al_Mytilus_spp_data.zip
Particle tracking.zip

EXPLANATION PER ZIP FILE
Coolen_et_al_Mytilus_spp_data.zip contains the following:
Coolen_et_al_Mytilus_spp_genotyping_data.csv; Which contains the genotyping results based on microsatellite data.
300_Ima2p_infiles; a folder containing all 300 Ima2p input files as csv and a single file containing list of file names as csv
300_Ima2p_outfiles; a folder containing all 300 Ima2p output files as csv

Particle tracking.zip contains the following:
2004 nieuw; a folder containing all particle tracking model results from the year 2004 as 25 csv files. File name relates to locations where particicles originate from. The order of particle receiving locations in the file (columns) is as follows:
HORN; BKUM; HARW; F31A; EURW; K10B; CALA; BG1x; PAAL; BLYT; EUR7; WZB; WZA; BRES; Q13A; G14A; K9A; L10G; L15A; BVW; D15A; FINO; SCHV; HELG; SYLT.
2005 nieuw; a folder containing all particle tracking model results from the year 2005 as 25 csv files. Otherwise identical to the 2004 nieuw folder above.

LOCATION INFORMATION
Location names and file names with location reference all refer to the following locations:
From left to right: Abbreviated name of sample, full name of sample; sampling date, structure type, sample size (number of specimens genotyped), sampling depth (m), position in decimal degrees WGS1984 (Lat. & Lon.) and distance to the nearest coastline (km). 
Name    Full name    Date    Structure    Size    Depth    Lat.    Lon.    Dist.
BG1x    BG 1    11 March 2015    Buoy    44    0    53.8833    3.4983    116
BKUM    Borkum Riffgat     29 June 2014    Wind farm    25    4    53.6900    6.4800    14
BLYT    Blythe    14 November 2014    Pier    48    0    55.1258    -1.4983    0
BRES    Breskens    20 August 2014    Breakwater    48    0    51.4068    3.5121    0
BVW    BV W    30 September 2014    Buoy    48    0    52.6007    3.5170    74
CALA    Calais    20 August 2014    Pier    24    0    50.9661    1.8433    0
D15A    D15-A    3 October 2015    O&G platform    48    7    54.3247    2.9346    181
EUR7    EURO 7    15 September 2014    Buoy    48    0    51.9900    3.5031    32
EURW    Euro W    24 July 2014    Buoy    48    0    51.9095    2.7232    70
F31A    F3-1A    1 September 2014    O&G platform    47    5    54.8520    4.6949    166
FINO    FINO 3    23 September 2015    Research platform    28    4    55.1950    7.1583    71
G14A    G14-A    16 July 2014    O&G platform    48    13    54.2241    5.4986    85
HARW    Harwich    15 November 2015    Pebble beach    42    0    51.9348    1.2813    0
HELG    Helgoland    15 January 2016    Harbour    48    0    54.1760    7.8945    47
HORN    Horns Rev    10 June 2015    Wind farm    67    0    55.4789    7.8110    19
K10B    K10-B    1 October 2014    O&G platform    48    27    53.3626    3.2539    104
K9A    K9-A    27 August 2014    O&G platform    48    0    53.5202    3.9925    66
L10G    L10-G    08 June 2014    O&G platform    48    10    53.4904    4.1952    53
L15A    L15-A    5 June 2014    O&G platform    48    6    53.3295    4.8302    11
LIMF    Limfjorden    16 June 2014    Longlines    48    2    56.7830    8.9110    0
LISB    Lisbon    14 February 2015    Harbour    37    0    38.7635    -9.0926    0
PAAL    Texel    29 June 2014    Breakwater    48    0    53.0118    4.7083    0
Q13A    Q13-A    28 May 2014    O&G platform    48    0 to 7    52.1911    4.1361    13
SCHV    Scheveningen    8 July 2014    Breakwater    48    0    52.0987    4.2582    0
SYLT    Sylt    3 June 2014    Breakwater    48    0    55.0216    8.4403    0
WZA    Wadden Sea A    30 April 2014    Intertidal mussel bed    24    0    53.4521    6.3042    0
WZB    Wadden Sea B    6 May 2014    Subtidal mussel bed    24    2    53.4600    6.3583    0