Marine stepping-stones: Connectivity of Mytilus edulis populations between offshore energy installations

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.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 Funding provided by: Dutch Department of Economic AffairsCrossref Funder Registry ID: Award Number: KB-24-002-001; KB-14-007Funding provided by: NWO Domain Applied and Engineering SciencesCrossref Funder Registry ID: Award Number: 14494Funding provided by: INSITE Foundation Phase ICrossref Funder Registry ID: Award Number: RECON & UNDINEFunding provided by: Nederlandse Aardolie Maatschappij BVCrossref Funder Registry ID: Award Number: Funding provided by: Wintershall Holding GmbHCrossref Funder Registry ID: Award Number: Funding provided by: Energiebeheer Nederland B.V.Crossref Funder Registry ID: Award Number: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 FST values directly estimated from the microsatellite data. Pairwise FST values were converted to 2Nem values following the infinite island model approximation by Wright (Wright, 1943) in which FST = (1 + 4Nem)-1. Correspondence between 2Nem values stemming from FST 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