Agent-based modeling and genetics reveal the Limfjorden, Denmark, as a well-connected system for mussel larvae

a well-connected ABSTRACT: Fishery of blue mussels Mytilus edulis constitutes a very important economic activity in Denmark, whereas mussel farming on long-lines or nets is a new, growing sector. Spawning from natural mussel beds takes place during early summer, and larvae disperse via water currents before settling on the bottom or on spat collectors in the farms. In the present study, we coupled a 3D physical model system (FlexSem) with an agent-based model in order to examine the connectivity of this marine system in terms of mussel larval dispersal and settling potential. To address this question, we (1) estimated the dispersal and connectivity between 17 areas in the Limfjorden, (2) identified the main donor and receiver areas of mussel larvae and (3) identified possible dispersal barriers. The results show that the central narrow strait in the Limfjorden was the main donor area in all the studied years, and that the adjacent eastern areas were the main receiver areas. Towards the inner basins of the Limfjorden, isolation increased and limited connectivity was observed. The results from the cluster analysis grouped the Limfjorden into 3 to 5 clusters, but there was still some exchange of simulated larvae observed among these clusters. Analysis of molecular markers revealed no genetic differentiation between areas and supports the model results, indicating that despite distinguishable hydrographic boundaries, the mussel populations in the Limfjorden are well connected. This study demonstrates how connectivity modeling can be used to support site selection processes in aquaculture.


INTRODUCTION
The marine environment provides different means for organisms to travel and disperse among populations. The life cycle of most marine organisms is characterized by a larval stage that is transported in the sea, leading to the potential of dispersal over broad geographic regions (Pineda et al. 2007, Cowen & Sponaugle 2009, Sala et al. 2013). This pelagic phase, which can last from a few days to several months, plays an important role in the connectivity of marine populations (Paris et al. 2007, Pastor et al. 2018. Connectivity is defined as the exchange of individuals among geographically distributed populations. Larval dispersal and population connectivity are highly linked because the processes that control dispersal of individuals from one location to another also connect marine populations demographically (Cowen & Sponaugle 2009). Quantifying the processes of dispersal and connectivity requires some degree of sim-plification due to the often complex and numerous physical and biological parameters involved on different spatial and temporal scales (Treml et al. 2008). Connectivity between 2 populations is dependent on the larval characteristics of the study species (e.g. pelagic larval duration, spawning month, vertical swimming behavior), the abundance of the source population, environmental factors (e.g. currents, temperature, salinity) and the availability of suitable habitat (Treml et al. 2008).
Invertebrates with pelagic larvae must successfully pass through several life cycle stages in order to ensure high recruitment to the population. Blue mussels Mytilus edulis produce large numbers of pelagic larvae that spend several weeks in the surface waters (Widdows 1991). Mussels begin spawning in early spring; this usually starts in March and peaks in April/May, depending on the water temperature (Coolen et al. 2020). M. edulis has a pelagic larval stage of 16−70 d, after which the larvae metamorphose to the pediveliger stage during which they are capable of settling (Coolen et al. 2020). The pediveliger larvae seek out different types of substrates before settlement, which involves some swimming and crawling behavior, culminating in the attachment of the larvae once a suitable substrate is chosen (Riisgård et al. 2015).
As a result of these complex cycles, many assumptions are made when studying larval dispersal with regard to larval supply, larval exchange, and larval spread. In order to overcome these assumptions, research efforts have been focusing on modeling environments coupling larval dispersal parameters with hydrodynamic data (Cowen et al. 2006, Paris et al. 2007, Treml et al. 2008, Hinrichsen et al. 2018). Agent-based models (ABMs) describe the behavior and state of individuals ('agents') or groups of organisms and simulate their responses in a spatio-temporal environment (Hansen et al. 2015). Agents are defined with a particular set of x-, y-and z-coordinates and behavioral mechanisms that can be linked to a hydrodynamic model (Hansen et al. 2015). This type of ABM is often referred to as a Lagrangian ABM and has been applied in different studies on marine invertebrates (Sale & Kritzer 2003, Yearsley & Sigwart 2011, Bendtsen & Hansen 2013 and fish larvae (Paris et al. 2007, Staaterman & Paris 2013, Daewel et al. 2015.
Developments over the last 25 yr in both physical oceanography and fish ecology have increased our ability to develop and integrate physical processes and biology (see review by Miller 2007). These concepts have contributed to the analysis of the dispersal of commercially valuable species such as cod in the Baltic andNorth Seas (Daewel et al. 2015, Huwer et al. 2016). Furthermore, they can be utilized to advise on the design of marine protected areas in different regions (Fogarty & Botsford 2007, Moksnes et al. 2014, Faillettaz et al. 2018.
Early studies on blue mussel larvae focused on physiological ecology (Sprung 1984, Widdows 1991, Van Haren & Kooijman 1993, first carried out in laboratories (Fotel et al. 1999, Pernet et al. 2003) and later during field campaigns (McQuaid & Phillips 2000, Dolmer & Stenalt 2010, Tilburg et al. 2012, Toupoint et al. 2012, to gain more information about larval settlement and seasonal variation in mussel abundance. In the following years, ABMs were used to investigate the dispersal of mussel larvae (Saraiva et al. 2014), with the aim of simulating population dynamics in the Wadden Sea. Genetic studies were also conducted to analyze Baltic populations of Mytilus and hybrid zones (Zbawicka et al. 2014). Recent studies have combined hydrodynamic modeling and genetics (Stuckas et al. 2017, Jahnke et al. 2020 in order to analyze the genetic composition of mussels in the North and Baltic Seas and to validate the models. The Limfjorden in Denmark is a large, shallow body of water lying across the Jutland peninsula with a number of shallow sub-basins and straits (Fig. 1). It is connected to the North Sea in the west and to the Kattegat (a strait between Demark and Sweden) to the east, but with limited water exchanges due to the narrow entrances (Hofmeister et al. 2009). There is a high standing stock of blue mussels in the area, and the mussel fishery is the largest in Denmark, with annual landings of 20 000−25 000 t wet weight (WW) for human consumption (data from the Danish Fisheries Agency, https://fiskeristyrelsen.dk/fiskeristatistik/ dynamiske-tabeller/). Suspended mussel aquaculture is a relatively new activity in the area, producing >2000 t WW for the fresh market (Danish Fisheries Agency), but is likely to expand in the future, since not all licenses are activated.
The optimal locations for mussel farms depend on local conditions such as current velocity, food concentration, salinity and temperature, but the supply of natural larvae to settle on the substrates in the farms is also important (Buer et al. 2020, Kotta et al. 2020, von Thenen et al. 2020. The supply of larvae to an area depends on the spawning events by the wild population, the transport of larvae by currents, the duration of the larval phase and the final settlement. However, no modeling studies have investigated how well the different basins of the Limfjorden are connected, whether dispersal barriers exist for mus-sels and which areas are the main suppliers for larval recruitment. Genetic analysis has shown that the mussels in the Limfjorden are primarily M. edulis (Kijewski et al. 2019). Further, 2 populations in the center and the north of the Limfjorden were similar to populations of the Kattegat, while a third population sampled in the south of the Limfjorden was similar to populations of the North Sea (Kijewski et al. 2019). Hence, it is possible that coexistence and hybridization of 2 differentiated genetic clusters could exist in the Limfjorden, possibly as a result of dispersal barriers and inflows from neighboring seas, and favor the detection of genetic structure and the analysis of connectivity with molecular markers (Gagnaire et al. 2015). However, no specific genetic studies have been conducted in this area.
In this study, we aimed to test whether the Limfjorden is a well-connected marine system for mussel larvae or if, in contrast, we can identify dispersal barriers limiting the transport of mussel larvae. To answer this question, we coupled the 3D physical model system FlexSem to an ABM in order to (1) estimate the blue mussel larval dispersal and connectivity between 17 areas (water bodies) within the Lim-fjorden, (2) identify main donor and receiver areas of mussel larvae and (3) identify dispersal barriers of mussel larvae in the system using both the model results and genetic studies. These findings can help us understand the dynamics in the system and can provide valuable information towards the optimal placement of mussel farms.

The Limfjorden regime
The total water area in the Limfjorden is 1575 km 2 , and the average water depth is 4.9 m. The center of the Limfjorden opens into the largest shallow basin called Løgstør Basin (Areas 10 and 11, see Fig. 1). Throughout most of the year, the wind blows from a westerly direction, with the exception of the summer period, which is dominated by easterly winds that are usually low in energy. The Limfjorden is a micro-tidal system with a tidal amplitude of 0.1−0.2 m, and it exchanges water with the North Sea (32−34 PSU, Area 1) and the Kattegat (19−25 PSU, Area 17) through narrow openings (Wiles et al. 2006). Furthermore, there is a freshwater input of 2.7 km 3 yr −1 from the catchment area. Currents, mixing and stratification in Limfjorden are governed by a complex interaction between wind-driven mixing and circulation, and density-driven currents (Hof meister et al. 2009).
The Limfjorden is a eutrophic water body affected by nutrient input from the surrounding watershed. This results in high primary production rates of up to 1000 mg C m −2 d −1 in summer (Maar et al. 2010), which supports a high level of biomass of benthic suspension feeders. The system is divided into 37 shellfish production areas ), but we have rearranged them into 17 areas for the purpose of the study (Fig. 1). A large proportion of these fisheries is located in Løgstør Basin (Areas 9−11). The majority of mussels are harvested from the bed, although mussels grown from long-lines are becoming more common (Taylor et al. 2019). There are presently 44 commercially licensed mussel farms in the Limfjorden, of which ~30 are currently active (Fig. 1). In recent years, cultivation practices designed Black dots represent the mussel bed sampling stations existing in the Limfjorden (used as release stations) and red dots indicate the current mussel farms. The system is divided into 17 areas for the connectivity study (area limits shown in blue). These areas have been modified from the current fishery areas for nutrient extraction have also been tested in the Limfjorden, with positive results (Nielsen et al. 2016, Taylor et al. 2019). There are currently 207 sampling stations for mussel stock as sessment ( Fig. 1), where mussel density (ind. m −2 ) has been recorded on a yearly basis since 1993 (Kristensen & Hoffmann 2004). The fishable area is covered by the sampling stations, consisting of squares of 0.7 × 0.7 nautical miles, and within each square dredging is conducted at a random position and direction at >3 m depth (Kristensen & Hoffmann 2004). Areas with stones and very few or no mussels, and areas protected from the fishery (some parts of Area 1, Area 4, the eastern part of Area 8, the northern part of Area 10, and Areas 16 and 17) were not sampled. We used the positions of the sampling stations as release areas in our model setup.

Hydrodynamic model
To model the physical environment in the Limfjorden, a 3D hydrodynamic model was made in the FlexSem framework (Larsen et al. 2020) using an unstructured computational mesh. The computational mesh combines triangles and squares to cover the Limfjorden by 6686 elements with an average horizontal resolution of 474 m (Fig. 2). The vertical resolution was in z-coordinates, i.e. the separation between computational cells in the vertical was defined at fixed depths, whereas the top layer had a free surface to allow for water level changes. The maximum water column depth was 30 m, resolved as 1.5 m in the surface layer, 1 m layer thickness in the upper 10 m, followed by 2 layers of 5 m each and a variable layer thickness at the bottom, to cover the remaining water depth. The surface layer was made sufficiently thick to accommodate the full range of surface level changes, thereby avoiding drying cells. Two open boundaries exchange water with the North Sea to the west and the Kattegat to the east (Fig. 2); hence, despite its name, the area is a sound, not a fjord. Initial fields and open boundary forcings of water level, horizontal velocities, temperature and salinity were obtained from nesting Limfjorden in the 3D ocean circulation Hiromb-BOOS model (HBM) (She et al. 2007, Berg & Poulsen 2012) run by the Danish Meteorological Institute. Likewise, 2D fields of atmospheric forcing of air temperature and wind speed and direction were interpolated from the HBM model output. River discharges of freshwater from 38 sources were obtained from the catchment model SWAT applied to the Limfjorden (Thodsen et al. 2016, Molina-Navarro et al. 2017. For more information on the hydrodynamic model validation see Supplement 1 at www.int-res. com/articles/suppl/m13559_supp.pdf.

ABM and parameter choice
To examine the connectivity in the Limfjorden and determine the potential for transport of mussel larvae, we used an ABM. In the ABM approach, individual agents are given an explicit x-, y-, z-coordinate at a given time and are transported by the currents. In each ABM time-step, and for each agent, the 3D Eulerian velocities are interpolated to the positions of the agents using an area-based interpolation approach, which then is advected in a Lagrangian way (Wang et al. 2011). For this study, an ABM module was coupled to the 3D hydrodynamic model from the Limfjorden setup (Fig. 2). Particle release was repeated for 5 yr (2009−2012 and 2017), representing years with a range of North Atlantic Oscillation (NAO) index values, which correlate well with the variability in circulation patterns (Kotta et al. 2020). The NAO index values were obtained from https://www.cpc. ncep.noaa.gov. More information about the ABM and its applications can be found in Larsen et al. (2020).
Mussel larvae in the ABM were defined by 2 biological parameters: the pelagic larval duration (PLD), which is defined as the period of larval development that is spent in the water column (Sponaugle et al. The PLD was set to 21 d, the reported value for mussel larvae (Widdows 1991, Riisgård et al. 2015, and the time of spawning was set to the beginning of May, which was found to be the peak month for mussel spawning in the Limfjorden by the Danish National Monitoring Database (https://odaforalle.au.dk). To simulate the dispersal of mussel larvae, numerical particles or agents, referred to here as 'simulated larvae,' were released from the mussel bed sampling stations acting as the source areas (Fig. 1). A constant number of larvae were released from all sampling stations (207 sites acting as larval sources), assuming an equal mussel density and spawning potential at those stations. As previously mentioned, some areas had no sampling stations, therefore, no larvae were released from those areas. A total of 80 000 simulated larvae were released during the main spawning event in May, which was subdivided into daily releases. The particles then moved freely in the water column after each release. The individual positions of the agents were updated every 3 min for the extent of the pelagic phase. A random speed of 0.01 m s −1 was incorporated in the agents' behavior (see Supplement 2) to account for water movement not resolved by the hydrodynamic model. Agents that hit land boundaries would be re-bounced to the center of the element in the mesh where they were in the previous time step. At the end of the PLD, the movement of larvae was deactivated, simulating larval settlement. The choice of parameterization is explained in the sensitivity an alysis in Supplement 2.

Connectivity analysis and clustering
Once the simulations were conducted, the connectivity of the system was calculated based on the downstream and upstream connectivity probabilities. Downstream connectivity is defined as the connectivity between a donor area and the different receiver areas. When simulating a large number of larvae, the equivalent large number of trajectories can be statistically analyzed, revealing the probability that an area will supply larvae to other areas (Hansen et al. 2015). This is referred to as a downstream connectivity probability. Upstream connectivity is de fined as the connectivity between a receiver area and the different donor areas. The probability that an area will receive larvae from other areas is referred to as the up stream connectivity probability. In this study, the areas refer to the 17 defined water bodies shown in Fig. 1. The connectivity probabilities (CP) are calculated as: ( 1) where X ij is the number of larvae going from area i to area j, and n is the number of areas. The downstream connectivity probability is therefore calculated as the number of larvae from donor area i that has settled in receiver area j divided by the total number of larvae originating from donor area i settled in all 17 areas n. Similarly, the upstream connectivity probability is calculated as the number of larvae in receiver area i that originate from donor area j, divided by the total number of larvae settled in receiver area i (Hansen et al. 2015). The resulting calculations are presented as connectivity matrices. Both downstream and up stream connectivity prob abilities include selfrecruitment values in the diagonal.
Following the connectivity study, we analyzed the topology of the transport network to subdivide the areas of the Limfjorden into clusters. The clusters and their boundaries are dynamical objects that evolve in space and time with different dimensions due to the important variability of the current circulation (Rossi et al. 2014). Based on the connectivity matrices, we detected separated clusters using the 'Infomap' algorithm (Rosvall & Bergstrom 2008). According to this method, agents are considered to move in a network system matching the statistical description of connectivity probabilities contained in the connectivity matrix. Using information theory concepts, Infomap decomposes the network into a number of communities that can define oceanic provinces (i.e. clusters of areas) well connected internally, but with minimal ex changes of larvae between them (Rossi et al. 2014). This clustering technique allows us to search for barriers in the dispersal. Little or no exchange be tween clusters could be an indication of possible dispersal barriers in the Limfjorden for the studied months.

Genetic analysis
The final stage of the study comprised the genetic investigation of mussel samples from different farms, in order to identify possible population differentiation in the Limfjorden and validate model results. The genetic study examined 23 single nucleotide polymorphisms al. 2020) and with a minor allele frequency of 0.05. The mussel samples were extracted from 7 farms in the Limfjorden, located in Venøsund, East of Venøsund, Sallingsund, Skive Fjord, Dråby Vig, Lovns Basin and Løgstør Basin (Fig. 1). In total, 40 mussels were processed per farm (except for Løgstør, where n = 26), and 10 mg of tissue were extracted from each mussel. Genotyping was subcontracted to LGC genomics and performed with the KASP™ array method (Semagn et al. 2014).
Once the genetic data were obtained, a principal component analysis (PCA) and a fixation index (F ST ) test for population differentiation were conducted in order to verify if the clusters from the connectivity analysis and the genetic samples matched. Individuals with more than 15% missing data were removed. Except for the PCA, mitochondrial markers were removed from the analyses. For the PCA, genotypes were centered and scaled using the R package 'adegenet' (v2.1.3, Jombart & Ahmed 2011), with the mean method for missing data replacement. We used the 'ade4' package (v1.7-15, Chessel et al. 2004) to compute the PCA. An analysis of molecular variance (AMOVA) (Excoffier et al. 1992) was computed using 10 4 permutations and the 3 levels of population, fishery area and cluster with the R package 'pegas' (v0.13, Paradis 2010). In this analysis, populations were grouped by predefined fishery areas and oceanographic clusters, based on their location regarding the ABM clustering results (see Figs. 1 & 5). F ST values (Weir & Cockerham 1984) and their confidence intervals (0.025 and 0.975 quantiles) were computed using a bootstrap with 10 4 permutations (R package 'hierfstat' v0. 04-22, Goudet 2005). Pairwise tests of significance using 10 4 permutations were also performed and followed by a Holm-Bonferroni correction for multiple testing. Finally, a modelbased genetic clustering algorithm was used to detect possible splits in the dataset. The program STRUCTURE (v2.3.4, Falush et al. 2003) was used and run with and without the admixture model. Twenty replicates were computed for each number of clusters (K) between 1 and 8, with a Markov chain Monte Carlo length of 80 000 steps after a 20 000-step burn-in. Results were aggregated with the program CLUMPAK (Kopelman et al. 2015) 3. RESULTS

Overall connectivity of the Limfjorden
Viewing the connectivity matrices allows one to discern spatial dispersal patterns. Fig. 3a,b shows the Fig. 3. Connectivity matrices indicating the (a) downstream and (b) upstream probability (proportion of larvae) for larvae originating in Areas 1−17 to end up in Areas 1−17 (indicated as 'donor' and 'receiver' areas). The diagonal elements indicate the probability of larvae staying in the same area where they were released (self-recruitment). Values are mean values from all studied years. Areas 4, 6, 16 and 17 had no (or only very few) mussel bed sites, and therefore no larvae were released from those areas downstream and upstream probabilities for all studied years. We can initially observe an overall high self-recruitment of simulated larvae in all of the areas in the Limfjorden across all years. This is represented in the diagonal values of the connectivity matrix (Fig. 3). Yet, apart from their own, mussel populations also seemed to be supported by other spawning areas nearby (Fig. 3b). Four out of the 17 areas located in the central and eastern part of the Limfjorden (4, 6, 16 and 17) did not have any or very few mussel bed sites, and therefore no larvae were released from those areas (Fig. 3a). Mussel larvae released from Areas 5 and 9, corresponding to the narrow straits in the central Limfjorden, represented donor larvae for most of the other areas in the west and north east (1, 2, 3, 6, 11 and 12) (Fig.  3a). These areas are also where most of the mussel farms are located. Area 13, which can also be defined as a strait, was further identified as a potential donor area to the surrounding water bodies, as evident in the downstream connectivity matrix (Fig. 3a). Areas 11, 12 and 13 in Løgstør received larval supply from many different areas (9, 10, 11 and 12). Very few farms are located in Løgstør (Areas 11 and 12). However, there are some in Area 13 (Fig. 1). We could also identify areas with few connections to the other water bodies, such as Lovns Basin and Skive Fjord (Areas 14 and 15), which were mainly supported by their own larval supply. Finally, the 2 eastern areas in Aggersund (Areas 16 and 17) did not provide larvae to other areas, since no mussels were released from them (Fig. 3a), although Area 16 could receive larvae from the neighboring Løgstør Basin (Fig. 3b) and is therefore not completely isolated from the system. In summary, the main donor areas in this study were primarily located in the narrow straits (Areas 5 and 9), the main re ceiver areas were in the central-eastern part (Areas 11−13), and the isolated areas were in the east and south-east of the Limfjorden (Areas 14 and 15 are primarily selfrecruiting and display little exchange, whilst Area 17 did not receive larvae from other areas).
The hydrodynamic data can help to explain the connectivity results obtained. Fig. 4 shows the mean and standard deviation of the current speeds in m s −1 within the simulation period (May−June) in the surface layer. The central areas of the Limfjorden, corresponding to Areas 5 and 9, show the highest changes in current speed and direction. Area 1 (Nissum Basin) also appears to have high mean current speeds. This is due to the tides that cause changes in current speed and direction within the mouth of the Limfjorden connected to the North Sea.

Cluster analysis
The Infomap algorithm used in this study identified 3−5 clusters, depending on the year analyzed (Fig. 5). The boundaries to these clusters seem to vary across years within the eastern areas and the central strait (Fig. 5). There was no larval release from Areas 16 and 17, so caution should be taken when differentiating those sites. The self-recruitment of larvae was >94% in 2009, 2010 and 2012 in all clusters, leaving little exchange to the rest of the clusters (Fig. 5). The years 2011 and 2017 presented lower self-recruitment rates of 73 and 74%, and therefore a higher exchange of larvae is observed among clusters for these years. Overall, clusters were rarely isolated from the rest, indicating that there was always an exchange of larvae and that the Limfjorden is well connected throughout the years.

Genetic results
We did not observe genetic clustering with the PCA (Fig. 6) or with the STRUCTURE model-based algorithm (Fig. S6 in Supplement 3). The AMOVA indicated that while the population level differentiation was significant (σ 2 = 0.26, df = 1, p = 0.03), the fishery area and cluster levels based on the oceanographic modeling results were not significant (p = 1 and 0.20, respectively). However, F ST values be tween all pairs of populations were non-significant after Holm-Bonferroni correction for multiple testing and all bootstrapped confidence intervals included 0 (Table S5 in Supplement 3). Observed and expected heterozygosities were similar between sites (Table S6 in Supple-ment 3). Overall, no genetic differentiation was observed between the populations, regions or oceanographic clusters of the Limfjorden.

Connectivity and mussel farms
The aim of the present study was to assess the mussel larval dispersal and connectivity in the Limfjorden, and to determine whether the system can be considered a well-connected system, or if, on the contrary, we can identify dispersal barriers limiting the transport of mussel larvae. Distinct dispersal barriers in the Limfjorden could have implications for the re cruitment of mussel larvae and optimal placement of mussel farms (Kotta et al. 2020). Overall, the results show that the self-recruitment of mussel larvae in all areas is generally high (>94% in 2009, 2010 and 2012, and between 70 and 80% in 2011 and 2017). However, mussel lar- vae are also dispersed to other areas, driven by the dominant currents during the spawning month.
The study also identified donor and receiver areas of mussel larvae (Fig. 3). The main donor area is the central strait in the Limfjorden corresponding to Areas 5 and 9 (Fig. 1), providing larvae to many other areas. This can be explained by the high mean and standard deviation of current speed observed for the period of the simulation (Fig. 4). This allows larvae to be transported in different directions and further into other areas before settlement. Interestingly, the straits are also where most mussel farms are located, suggesting that this area supports a high food flux for mussel farming (Taylor et al. 2019). The main receiver areas are located on the eastern side of the Limfjorden corresponding to Areas 11, 12 and 13 (Fig. 1), and again, this is where a high concentration of mussel farms is located. These mussel farms might benefit not only from the high self-recruitment, but also from a sustainable supply of larvae from other areas making them less vulnerable to local changes in recruitment. As a last remark regarding the overall connectivity, the 2 south-eastern inner areas of Lovns Basin and Skive Fjord (Areas 14 and 15) could be defined as the most isolated parts of the system that are primarily self-sustained. The absence of farms in this isolated area can be explained by the periods of severe hypoxia that occur every summer (Møhlenberg 1999), where mussels cannot survive in the stratified parts and therefore do not provide suitable conditions for mussel growth and recruitment.
This study primarily focuses on potential connectivity, as defined by Watson et al. (2010), as opposed to realized connectivity, which attempts to quantify the absolute numbers of larvae that connect different locations. Potential connectivity is a relative estimate used to identify which regions may be connected using relative terms. In the absence of data on mussel densities from sampling stations and of data on predicted mussel habitat distributions, the distribution of sampling stations is considered representative of the spatial distribution of mussels in the Limfjorden. Alternatively, if data had been available, the study could have been based on a more accurately predicted mussel distribution. However, given the high connectivity in the entire system and the high density of stations in the central parts of the Limfjorden, where fisheries and mussel cultures are present, this would presumably have led to similar conclusions.

Clustering and genetics
The current model study identifies 3−5 clusters in the Limfjorden when applying the Infomap algorithm (Fig. 5). These clusters are pooled areas, which can be differentiated from the rest and present high selfrecruitment. However, no cluster is completely iso-9 Fig. 6. PCA plot with the 2 oceanographic clusters predicted by the model (C1 and C2) and the populations (the areas where the mussel samples were taken, see Fig. 1). Principal components 1 and 2 (PC1 and PC2) explain 6.7 and 6.4% of the variation in the data, respectively lated from the other clusters as seen from the connectivity results (Figs. 3 & 5), and there is also some yearto-year variability in their geographical coverage, suggesting a well-connected system. The genetic results support this idea and did not reveal any differentiation between the analyzed sampling sites, as shown in the PCA. Contrary to the results of Kijewski et al. (2019), our samples from the south of the Limfjorden were not genetically differentiated, and our results support genetic panmixia at the scale of the whole study area. Note that the genetic composition of Mytilus edulis in the Limfjorden is more similar to those in the Kattegat and Skagerrak than to M. edulis from the North Sea, as was found by Kijewski et al. (2019), but only 1 genetic cluster was observed in our dataset. We suspect that a transition zone might exist between populations of the Limfjorden and the North Sea, west of the VE site (potentially in Area 1, see Fig. 1), and that Kijewski et al. (2019) sampled the North Sea lineage while we did not. Here, we address the genetic differentiation within the Limfjorden lineage which proved panmictic. A small level of migration in each generation is sufficient to effectively homogenize the genetic composition of large marine populations in the absence of processes that maintain disequilibrium (Gagnaire et al. 2015). In addition, mussels have a high reproductive output, and although there is high self-recruitment in the identified clusters, the currents in the system still allow larvae to reach neighboring areas or in some cases to be transported even further away. This exchange seems to be enough to support the connectivity among areas, and smooth out potential genetic differences resulting in a well-mixed genetic pool. We should note that mussels were collected in farms and this could potentially influence our results in 2 ways. First, while cultivation practice suggests that mussel spat is collected locally by mussel farmers, we cannot totally exclude that spat is not sometimes imported from another area. Second, a different genetic composition might exist between natural mussel beds and rope cultures due to settlement and post-settlement selection, especially at the entrance of the Limfjorden where 2 lineages might co-exist according to our results and those of Kijewski et al. (2019). Differential recruitment selection at small spatial scales is well known in zones where mussel lineages occur in sympatry (Comesaña & Sanjuan 1997, Wilhelm & Hilbish 1998, Katolikova et al. 2016. The results of genetic and connectivity studies can be considered to be in good agreement in this study, and indicate that self-recruitment of mussel larvae should be >94% before we observe clusters as actual dispersal barriers in the Limfjorden. The comparison between genetics and connectivity assessments has been conducted in other studies with very different outcomes. Jahnke et al. (2016) assessed long-distance dispersal of the seagrass Zostera noltei in the Black Sea and found good agreement between both methods. Good concordance was also observed in Z. marina in the Skagerrak−Kattegat region and the Swedish west coast, but only with multigenerational models (Jahnke et al. 2018(Jahnke et al. , 2020. In contrast, Johansson et al. (2015) investigated the population genetic structure of the giant kelp Macrocystis pyrifera in the northeast Pacific and found that the genetic clusters were not supported by oceanographic transport. In mussels, although Gilg & Hil bish (2003) managed to link genetics and hydro dynamic models in a hybrid zone maintained by a balance between migration and reproductive isolation mechanisms, the approach is likely to be much less successful within homogeneous genetic clusters at equilibrium. Recently, Coolen et al. (2020) investigated the marine stepping-stone effect on M. edulis in the North Sea, and found no genetic differences despite the distinct dispersal barriers found in the connectivity modeling study. Barriers to dispersal were probably as efficient in the North Sea offshore platforms as in the Cornwall hybrid zone, if not more, but were likely insufficient to maintain genetic differentiation in a context of genetic equilibrium (Gagnaire et al. 2015). Mussels could also be effectively dispersed by other vectors such as fishing vessels, leisure boats, or floating macro algae (especially in the Limfjorden) and in this way, hide the 'signal' of hydrographic barriers in the genetic analysis.
Overall, these studies with different outcomes highlight the complexity and importance of using both methods in order to understand species connectivity in different marine systems.

Implications for marine management
In this study, we show that modeling tools combined with species genetics can provide detailed information on the dispersal and potential connectivity of a species in an area. Both approaches showed that Limfjorden is a well-connected system for mussel larvae dispersal and recruitment. Questions towards finding suitable places to establish a mussel farm might be challenging, as many aspects like mussel growth, food supply, predation and the risk of low oxygen conditions have to be considered (Friedland et al. 2019). However, the present findings can help us understand the dynamics of the system and support the identification of areas with high potential for mussel release and settling, which could be utilized for future mussel farming.
Gaining better understanding of the mussel recruitment processes can help us adapt to climate change and better solve environmental management questions in different areas. It also raises the possibility of identifying areas where protection of mussel populations is important, and which areas could be fished more intensively with limited impact on larval recruitment and connectivity in the system. Habitat suitability models have been used to identify suitable locations for mussel farming from a spatial and environmental perspective (von Thenen et al. 2020). These models can also be combined with ABMs in order to support management decisions and marine spatial planning.
Further work will include investigation of changes in climate, which are expected to affect the larval dispersal probabilities in the system through the currents and the biological characteristics of the mussels. Our framework has direct applicability to many other marine systems with aquaculture, and may assist in gaining a more holistic and integrated view of dispersal-based connectivity to aid management in the site-selection process with respect to recruitment of larvae and optimal placement of shellfish farms.