DOI: 10.5061/dryad.612jm6405
title: Marine stepping-stones: Connectivity of Mytilus edulis populations between offshore energy installations

description: 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 &gt;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.
description: <i>Mytilus
 edulis</i> 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
 <i>M. edulis</i> (Denmark) and <i>M.
 galloprovincialis</i> (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
 <i>M. edulis</i>, <i>M. galloprovincialis</i> and
 <i>M.</i> <i>trossulus</i> (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 (<i>e.g.</i>,
 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 <i>M. edulis</i>.
 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 <i>M. edulis
 </i>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 <i>Mytilus galloprovincialis</i> and <i>M.
 trossulus</i> 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 (<i>M.
 galloprovincialis</i>) and Limfjorden (<i>M. edulis)</i>
 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 <i>M. edulis</i>
 (<i>h</i> &lt; 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 <i>t</i>, population sizes
 <i>q</i> and bidirectional migration rates
 <i>m</i> 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
 <i>q</i> (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 <i>t</i> (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<sub>ST</sub> values directly
 estimated from the microsatellite data. Pairwise
 F<sub>ST</sub> values were converted to
 2N<sub>e</sub>m values following the infinite island model
 approximation by Wright (Wright, 1943) in which F<sub>ST</sub>
 = (1 + 4N<sub>e</sub>m)<sup>-1</sup>.
 Correspondence between 2N<sub>e</sub>m values stemming from
 F<sub>ST</sub> 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.
description: EXPLANATION OF DATASET<br> The data are formatted in csv
 files contained in two zip files:<br>
 Coolen_et_al_Mytilus_spp_data.zip<br> Particle
 tracking.zip EXPLANATION PER ZIP FILE<br>
 Coolen_et_al_Mytilus_spp_data.zip contains the following:<br>
 Coolen_et_al_Mytilus_spp_genotyping_data.csv; Which contains the
 genotyping results based on microsatellite data.<br>
 300_Ima2p_infiles; a folder containing all 300 Ima2p input files as csv
 and a single file containing list of file names as csv<br>
 300_Ima2p_outfiles; a folder containing all 300 Ima2p output files as
 csv Particle tracking.zip contains the
 following:<br> 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:<br> 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.<br> 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<br> Location names and file names with location
 reference all refer to the following locations:<br> 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). <br> Name    Full name  
  Date    Structure    Size    Depth    Lat.    Lon.    Dist.<br>
 BG1x    BG 1    11 March 2015    Buoy    44    0    53.8833    3.4983  
  116<br> BKUM    Borkum Riffgat     29 June 2014    Wind farm  
  25    4    53.6900    6.4800    14<br> BLYT    Blythe    14
 November 2014    Pier    48    0    55.1258    -1.4983    0<br>
 BRES    Breskens    20 August 2014    Breakwater    48    0    51.4068  
  3.5121    0<br> BVW    BV W    30 September 2014    Buoy    48  
  0    52.6007    3.5170    74<br> CALA    Calais    20 August 2014  
  Pier    24    0    50.9661    1.8433    0<br> D15A    D15-A    3
 October 2015    O&G platform    48    7    54.3247    2.9346  
  181<br> EUR7    EURO 7    15 September 2014    Buoy    48    0  
  51.9900    3.5031    32<br> EURW    Euro W    24 July 2014  
  Buoy    48    0    51.9095    2.7232    70<br> F31A    F3-1A    1
 September 2014    O&G platform    47    5    54.8520    4.6949  
  166<br> FINO    FINO 3    23 September 2015    Research platform  
  28    4    55.1950    7.1583    71<br> G14A    G14-A    16 July
 2014    O&G platform    48    13    54.2241    5.4986  
  85<br> HARW    Harwich    15 November 2015    Pebble beach    42  
  0    51.9348    1.2813    0<br> HELG    Helgoland    15 January
 2016    Harbour    48    0    54.1760    7.8945    47<br> HORN  
  Horns Rev    10 June 2015    Wind farm    67    0    55.4789    7.8110  
  19<br> K10B    K10-B    1 October 2014    O&G platform  
  48    27    53.3626    3.2539    104<br> K9A    K9-A    27 August
 2014    O&G platform    48    0    53.5202    3.9925  
  66<br> L10G    L10-G    08 June 2014    O&G platform  
  48    10    53.4904    4.1952    53<br> L15A    L15-A    5 June
 2014    O&G platform    48    6    53.3295    4.8302  
  11<br> LIMF    Limfjorden    16 June 2014    Longlines    48    2  
  56.7830    8.9110    0<br> LISB    Lisbon    14 February 2015  
  Harbour    37    0    38.7635    -9.0926    0<br> PAAL    Texel  
  29 June 2014    Breakwater    48    0    53.0118    4.7083    0<br>
 Q13A    Q13-A    28 May 2014    O&G platform    48    0 to 7  
  52.1911    4.1361    13<br> SCHV    Scheveningen    8 July 2014  
  Breakwater    48    0    52.0987    4.2582    0<br> SYLT    Sylt  
  3 June 2014    Breakwater    48    0    55.0216    8.4403    0<br>
 WZA    Wadden Sea A    30 April 2014    Intertidal mussel bed    24    0  
  53.4521    6.3042    0<br> WZB    Wadden Sea B    6 May 2014  
  Subtidal mussel bed    24    2    53.4600    6.3583    0