Integrating seascape resistances and gene flow to produce area-based metrics of functional connectivity for marine conservation planning

Prioritizing regions that facilitate connectivity among populations is an essential principle for conservation planning. However, the lack of conspicuous geographical and environmental features that constrain dispersal and gene flow throughout life history challenges the characterization of dispersal pathways within a three-dimensional marine realm. To elucidate regions of high connectivity value in the marine environment, we develop a novel approach that integrates estimates of spatial genetic structure with representation of regions of high dispersal potential for meroplankton, incorporating elements of pelagic larval and benthic adult life history. Spatial patterns of connectivity were characterized using circuit theory as an inverse function oceanographic- and habitat-based resistance to movement. We integrate emergent spatial patterns of connectivity with population genetic data to account for realized patterns of gene flow across a seascape. We apply this approach to four broadly distributed species in the Northwest Atlantic. Estimates of resistance to gene flow revealed multiple connectivity barriers not observed in oceanographic or habitat models. Comparison of isolation-by-distance versus isolation-by-resistance revealed genetic variation was best explained by seascape resistance in three of four species, supporting the resistance-based assessments of connectivity. Our approach identified areas of high and low connectivity value for each species, with overlap generally associated with geographic pinch points and areas of low genetic exchange. By integrating spatial interpolations of gene flow and estimated pathways for dispersal, we develop a novel area-based metric of connectivity that considers life-history based structural constraints to dispersal and observed genetic variation. Outputs from this workflow can reveal regions of connectivity for conservation planning.


Introduction
Connectivity is a fundamental ecological process that involves the flow of nutrients, organisms, genes, and/ or energy among spatially and/or temporally distinct entities (Kool et al. 2013). Connectivity strongly influences population, community, and ecosystem structure ) and can promote persistence and recovery in the face of natural and anthropogenic stressors (Burgess et al. 2014). From an organismal and/or genetic perspective, this broadly defined process is commonly partitioned into demographic, genetic, structural (herein seascape), and functional connectivity. Seascape connectivity concerns how physical characteristics of the environment (e.g., temperature, salinity, circulation, and seafloor topography) impede or promote movement of organisms (Taylor et al. 1993). Demographic or population connectivity concerns the exchange of adults, larvae, and/or propagules among metapopulations, which influences demographic rates. Genetic connectivity concerns the dispersal of genes, which unlike demographic connectivity, accounts only for dispersive individuals that contribute to gene flow (i.e., those that move between locations and successfully reproduce). Finally, functional connectivity considers the outcome of species' interactions with the environment during recruitment and dispersal, reflecting seascape, demographic, and genetic connectivity (e.g., Metaxas and Saunders 2009).
Connectivity has long been a central component of conservation planning in terrestrial systems (e.g., Razgour et al. 2014) and recently has begun to be applied in the marine setting. With growth in the global protected area network driven by international conservation commitments (i.e., Aichi Target 11), there is a growing community of science working to understand how the design of protected areas, in particular protected area networks, can facilitate persistence of species and ecosystems (e.g., Magris et al. 2016;Balbar and Metaxas 2019;Friesen et al. 2019). Measures to protect and restore marine habitats such as the global conservation targets, explicitly acknowledge both the need to conserve diversity through ecological representation but also through establishing well-connected systems of protected areas (CBD 2010). The incorporation of connectivity into marine spatial planning and conservation design can serve to influence the dynamics of the species and environments represented and to promote their resilience (Burgess et al. 2014). From a demographic perspective, connectivity among protected areas can be particularly important when the spatial scale(s) of dispersal precludes single areas facilitating selfrecruitment and independently viable populations (Marti-Puig et al. 2013). Similarly, from a genetic perspective, connectivity among metapopulations can facilitate an increased capacity for genetic adaptation in a changing environment (Xuereb et al. 2019).
Though connectivity has often been identified as a key component of marine protected area (MPA) planning (e.g., Balbar et al. 2020), the realized application of connectivity in the design of MPA networks remains an area needing improvement (Magris et al. 2014;Balbar and Metaxas 2019). An MPA network is a collection of spatially discrete MPAs that promote and provide protection for ecological connectivity processes among the constituent areas. Though the development of MPA networks has increased globally, few meet the definition of a connected network (Grorud-Colvert et al. 2014). Currently, there are minimal practical applications of the research from scientific studies that evaluate connectivity into MPA network design, with connectivity being among the least frequently used ecological criteria (11% of MPAs examined-Balbar and Metaxas 2019), despite its perceived importance. This lack of incorporation could serve to jeopardize the long-term viability and effectiveness of MPA networks globally (Magris et al. 2018), particularly as ecosystems shift and respond to global climate change Wilson et al. 2020;Bryndum-Buchholz et al. 2022).
There remain challenges with incorporating connectivity into marine conservation planning. With such a broad definition, it is unsurprising that there are myriad metrics and methods, relating to the various types of connectivity, available for both scientists and practitioners to evaluate connectivity patterns (Botsford et al. 2009;Bryan-Brown et al. 2017;Balbar et al. 2020). The ineffective implementation of connectivity in MPA network design globally could, in part, stem from a disconnect between scientists and practitioners working on MPAs, with respect to some disparity and bias for examining certain types of connectivity (Balbar and Metaxas 2019). Ultimately, the appropriateness of any measure will depend on available data and the established conservation priorities (Smith and Metaxas 2018); however, robust assessments of connectivity should incorporate multiple metrics that examine the various types of connectivity at ecologically relevant spatial scales. Approaches that include complimentary measures can provide a more comprehensive understanding of connectivity between current and prospective conservation areas.
The area-based focus of traditional conservation planning tools may inadvertently contribute to the lack of incorporation into network design, where many metrics of connectivity are not represented spatially across the planning landscape. Various metrics of geographic and genetic distances have been used to identify important environmental variables or physical distances that structure populations and constrain dispersal (Balbar et al. 2020;Bertola et al. 2020). However, these metrics generally lack information detailing the particular value of areas to connectivity. Network metrics such as centrality and local retention (e.g., D'Aloia et al. 2017) offer another approach for representing and incorporating connectivity into decision making . Habitat suitability models provide area-based metrics that can be implemented to evaluate (relative) transmittance (probability of movement based on habitat preference) and simulate how dispersal and connectivity would be distributed across a landscape (Razgour et al. 2014;Barbosa et al. 2018). Biophysical-and habitat-based approaches are informative of potential connectivity, indicating where dispersion is likely to occur and/or be successful. For marine species with dispersive larvae, genetics can ground truth dispersal models and reveal important information on post-settlement processes, which often are an important driver of demography and distribution but are otherwise missed by studies which only characterize pelagic connectivity (Mertens et al. 2018). This genetic structure can be interpolated spatially. For example, Estimated Effective Migration Surfaces (EEMS) represent gene flow over a defined area as a function of the pairwise genetic dissimilarity among sampled sites (Petkova et al. 2016). When matched to environmental gradients (like in our study region- Stanley et al. 2018), the EEMS output is particularly advantageous as it provides a parsimonious, spatial-representation of gene flow representative of key environmental barriers to connectivity. Combining spatial habitat-based resistance, oceanographic dispersal, and gene-flow would provide a novel, comprehensive, and explicitly spatial representation of functional connectivity that could be applied to conventional area-based marine conservation planning and directly incorporated into planning software like Marxan.
In this study, we propose, develop, and test a methodological framework for assessing functional connectivity in a marine environment through an integration of seascape and genetic connectivity. To accomplish this, we develop approaches that integrate spatial genetic variability and habitat resistance to identify regions of connectivity through a planning seascape, building on approaches that to date have been exclusive to the terrestrial realm (e.g., Barbosa et al. 2018). We apply this approach using four ecologically and commercially significant marine species from the Northwest Atlantic, collectively representing a variety of life history strategies, habitat preferences, and scales of larval and adult movement, providing a broad representation of dispersal phenotypes for benthic fish and invertebrates present in the study area.

Materials and methods
The objective of our approach was to integrate information on habitat, oceanographic and genetic resistance to dispersal to illustrate regions important for connectivity across the planning area. Our study area encompasses northeastern United States and Atlantic Canada (Fig. 1). Steps for our analysis are described sequentially (Fig. 2).

Population genetic structure
We used published genetic data (Table 1) for Atlantic Cod (Gadus morhua), American Lobster (Homarus americanus), Atlantic Sea Scallop (Placopecten magellanicus), and European Green Crab (Carcinus maenas) to estimate gene flow among sampling sites and compare relationships of isolation-by-distance (IBD) against isolation-by-resistance (IBR). We used the R package genepopedit  to select the putative, neutrally evolving SNPs by subsetting outlier SNP panels identified in previous studies from full SNP datasets. Arlequin v.3.5.2.2 (Excoffier and Lischer 2010) was used to calculate pairwise linearized F ST (Slatkin and Voelm 1991) among populations from the putatively neutral SNPs. When genotyping a sufficient number of genome-wide SNPs (i.e., > 1000), F ST can provide a robust metric for characterizing genetic differentiation of marine populations, even when sample sizes themselves are small (Willing et al. 2012). Only neutral SNPs were used to calculate F ST and estimate gene flow to remove the influence of potentially adaptive (outlier) SNPs that can skew migration estimates, and because the inclusion of outlier SNPs resulted in a non-Euclidean dissimilarity matrix for estimating migration surfaces in EEMS (i.e., the eigenvalues for the matrix sum did not sum to zero as required in a Euclidean matrix; Petkova et al. 2016).
Gene flow across the seascape Effective migration surfaces were used to evaluate geographic regions of higher-or lower-than-average gene flow for each species. This method requires a user-defined number of demes to approximate all the possible paths between two individuals based on genetic ancestry, using resistance distances based on circuit theory (Petkova et al. 2016). For this study, 1000 demes provided a suitable scaffold (allowing passage and connection of demes through narrow straits and bays while preventing passage across land; Online Resource 1) to interpolate migration for each Fig. 1 A map of the Northwest Atlantic study area with relevant area labels. The range of genetic and habitat suitability data examined for the species in this study generally ranges from the northeastern United States to Newfoundland. The inset shows the broader North Atlantic Ocean and landmasses species. Missing loci were imputed as the average genotype at that SNP and dissimilarity matrices were created from the genotype files using the str2diffs function (https:// github. com/ dipet kov/ eems/). The resulting genetic dissimilarity matrices, a coordinates file for sampling locations, and a geographic bounding polygon where land was presented as a barrier to dispersal, were used as input in the software EEMS (Petkova et al. 2016). EEMS was run in triplicate for each species, with a burn-in of 200,000 iterations, Fig. 2 Flow diagram depicting the steps building towards a geographic representation of functional connectivity. Example rasters based on Atlantic scallop are meant to be emblematic of the integration of spatial information. Dashed lines represent the combination of two layers, with arrows pointing to the product. For functional connectivity, EEMS layers are multiplied by the optimized Circuitscape resistance chosen by the maximum correlation between pairwise Circuitscape resistance (CC res ) and linearized F ST for each species (Table 2). Each map depicts the study domain and the districts used to calculate current bearing are noted in the upper left panel. Legends depict colour scales of high to low resistance to movement, whereby high functional connectivity occurs in areas of low resistance Markov Chain Monte Carlo (MCMC) length of 2,000,000 iterations, and a thinning every 1000 iterations. Maps and associated rasters of effective migration rates were generated (Reemsplots, R), with rasters subsequently used when calculating functional connectivity in combination with seascape connectivity output(s) (see below- Fig. 2). While EEMS provides estimates of gene flow across our study area, it does not explicitly take habitat features or oceanography into account (aside from defined barriers such as land), which are needed to characterize functional connectivity.

Circuitscape resistances
We examined patterns of IBD and IBR by comparing metrics of geographic distance and landscape resistance (respectively) to genetic differences (F ST ) between sites for each species. To evaluate an initial approximation of the IBD relationship, we calculated pairwise great circle distance between sample sites in R. To examine relations of IBR, we developed resistance layers relative to adult (based on habitat suitability modelling) and larval (based on a simple advective current vector layer) dispersal. Given that dispersal over successive generations is a product of both adult movement and larval dispersal, we developed an additional layer that integrated both seascape and oceanographic resistance, representing habitat constraints for connectivity over lifespans. This layer represents both the potential dispersal of planktonic larvae and adult habitat necessary to support intermediate populations that would facilitate connection between our sampling sites over successive generations. In each surface, pairwise resistance among populations was evaluated using Circuitscape v5 (McRae 2006), which relates effective resistance 'distances' among sites (nodes through which electrical current travels between) precisely to random walks, such that circuit theory is related to movement ecology across complex landscapes (McRae et al. 2008). Compared with conventional methods that use least-cost paths or straight-line, Euclidean distances to assess spatial relations with gene flow (i.e., IBD), the use of resistance surfaces in Circuitscape assumes species are capable of using the entire seascape allowing for the assessment of dispersal over multiple pathways, which not only produces a metric for connectivity but can be used for identifying broad regions of high connectivity and important features of the seascape over which dispersal occurs (McRae 2006;McRae et al. 2008). Additionally, when trying to incorporate larval dispersal processes, the uncorrelated and unbiased transmission of current, which approximates a simple random walk model, relates to diffusive processes (Weiss 1994;Codling et al. 2008), while the resistance layer of ocean circulation incorporates advective processes. Using the 'pairwise mode' in Circuitscape, effective resistance values and conductance maps were generated for each site combination. These conductance maps depict the probability of movement through the seascape between site pairings, which were combined into a cumulative map to visualize movement between all sites within the study area.

Habitat resistance
To evaluate habitat-based resistance, we employed previously published population centric habitat suitability predictions for our study domain that were developed to account for environmentally associated population structure (northern and southern 'ecotypes') at a resolution of ~ 1 km . Suitability model predictions for the study area were transformed into resistance layers by rescaling to a range between 0.0001 (high suitability, low resistance) to 1 (low suitability, high resistance) according to McRae et al. (2008). Values for unsuitable marine habitat were set at 2, representing the increased difficulty of moving through or inhabiting a particular area. Terrestrial areas within the study domain were masked and thus treated as impermeable barriers for transmittance. Conductivity (i.e., the inverse of resistance) analyses were run for each resistance layer for each species ecotype. For each site pairing, the lowest effective resistance value between the ecotype resistances generated by Circuitscape, was selected.

Oceanographic resistance
Resistance layers of ocean circulation were developed using similar methods to Dambach et al. (2016), where current direction represents a potential resistance to connectivity (i.e., connection between sites likely to be where the current direction facilitates dispersal). Values for resistance were calculated based on the deviation of current direction from the bearing angle to the destination site. A deviation value of 0° represents no resistance (current direction = bearing angle) and a deviation value of 180° represents maximum resistance, which would be standardized for resistance through division of the maximum potential deviation (e.g., 180°). The output from this resistance layer analysis represents a simplified model of advective dispersal based on current bearing. Resistance layers were produced for each site considering the bearing from each pairwise site (acting as the supply of current) to the focal node (the site for which electrical current is drawn to) and thus are both unique for each focal node and asymmetrical within pairings (i.e., the resistance from site A to B ≠ resistance from site B to A). To develop the oceanographic resistance layers, we extracted surface current vectors (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) from the Bedford Institute of Oceanography North Atlantic Model (BNAM- Wang et al. 2018). For each species, currents were averaged for periods that encompassed expected spawning seasons and periods of expected larval pelagic dispersal (Online Resource 2). Average vectors were used to determine the deviation from bearings to the focal node. Unlike in Dambach et al. (2016), whose models of resistance were established for open ocean, the coastal environment of our study contains impassable landforms that sometimes negate the use of direct bearings. To account for these landforms, we divided the study area into "districts" where deviations from bearings were based on waypoints that lead to the boundary of the adjacent district closest to the focal node (Fig. 2). This method facilitated navigation around impassible landforms between sites. In one district located at the entrance to the Gulf of Saint Lawrence, the bearing direction for certain focal node sites lead to two possible directions to travel around the island of Newfoundland (Fig. 1). In such cases, the resistance values for this district were calculated using a median bearing direction between the two desired bearings (e.g., if bearings were 270 and 90, then bearing was set at 0). Resistance was then calculated as the absolute difference in the current deviation from the deviation of the median bearing to the desired bearings and standardized by dividing by the greatest potential deviance (e.g., resistance =|deviation-90|/90). Each resistance layer was produced for a grid at a nominal resolution of 7.5 km, where the values of each cell represented the mean resistance for all stations within each cell. Data was further disaggregated by a factor of 8 to facilitate the masking of land barriers, while ensuring connection through narrow passages of water. The known asymmetrical nature of the current resistances based on focal location resulted in two values for each pairing. Matrices of effective resistance were therefore constructed as an aggregate of values using the shortest of these two distances for each pairing of sites.

Composite resistance
To reflect how oceanographic currents and habitat act to jointly constrain dispersal and potentially gene flow across planktonic and benthic life stages, we developed a composite resistance layer that integrates habitat-based and oceanographic resistance layers. Composite resistance layers were produced as the product of habitat (each ecotype) and prevailing current resistance (each focal site) for each species at a resolution of 1 km. Resistance values of 0 were reassigned to a value of 0.0001 and a land mask was similarly applied to each raster prior to running through Circuitscape. As with the previous examinations, construction of resistance matrices considered only values derived from the specific ecotype when pairings of sites were contained within that ecotype (Fig. 2). Considering the asymmetrical nature of the oceanographic rasters, resistance values were selected based on the lowest effective resistance value between site pairings.

Functional connectivity
Seascape representations of functional connectivity were generated as the product of resistance and gene flow. First, to evaluate which resistance metric (habitat, oceanographic, or composite) best reflected gene flow or connectivity across the study area, we compared relationships between pairwise effective resistance metrics (and geographic distance) and population-level genetic differences (linearized F ST ) using Pearson's product-moment correlation. Cumulative conductance maps for the resistance metric that best reflected gene flow were generated through the addition of effective resistance values for each specific pairwise comparison used to assemble aggregate matrices of those metrics. Cumulative conductance maps were then rescaled to values between 0 and 1, with values representing habitat with low and high connectivity, respectively. Surfaces from EEMS were also standardized to 0 and 1, reflecting lower and higher than average migration. Final maps of functional connectivity for each species were generated as the product of the cumulative current maps and the effective migration surfaces using ArcMap 10.4.

Seascape genetics
All relations of IBD and IBR revealed significantly positive correlations between connectivity metrics and genetic differentiation (F ST ) except for the comparison of resistance based on prevailing current for lobster populations (p = 0.152) ( Table 2, Online Resource 3). For cod, the highest correlation resulted from comparisons with effective resistance from habitat resistance layers (r = 0.828). Lobster exhibited the highest correlation when comparing genetic differentiation with great circle distance (r = 0.418); however, for the purpose of evaluating an area-based metric, habitat resistance layers (r = 0.335), the next highest correlation, were selected as the metric to produce functional connectivity maps. Both crab and scallop exhibited the highest correlation when comparing the effective resistance from composite resistance layers (r = 0.704 and 0.588, respectively).

Estimates of gene flow
For all four species, the EEMS results agreed with the previously reported population structure, separating populations into northern and southern ecotypes with a region of low migration between them associated with a steep thermal gradient ) (Fig. 3). In addition to the primary region of reduced gene flow, Atlantic cod along the northern Scotian Shelf exhibited a potential gene flow barrier between populations from the Gulf of Saint Lawrence and eastern Newfoundland. Lobster exhibited the weakest barriers among populations, except for the Boothbay population in the Gulf of Maine, which showed reduced migration rates with, and relatively high isolation from, all populations in its vicinity. In sea scallop, Little Bay in southern Newfoundland showed reduced migration to other northern populations, as reflected by its observed population structure (Van Wyngaarden et al. 2017) (Fig. 3). In green crab, the Kejimkujik National Park population in southern Nova Scotia is a hybrid population between northern and southern ecotypes (Jeffery et al. 2017b). Since this hybrid population contains an approximately equal number of northern and southern ecotype alleles, it shows strong barriers to gene flow immediately to the north and south, likely reflective of the genetic break between the "pure" northern and southern ecotypes. Among populations within each northern and southern ecotype division, gene flow was relatively high, as would be expected (Fig. 3).

Seascape resistance maps
Cumulative conductance maps were assembled based on the highest significant correlation among the resistance metrics tested for each species, with cod and lobster based on habitat resistance and crab and scallop based on a composite resistance. Areas important for connectivity (i.e., areas with high current density) of the more mobile and wider ranging species tended to be more dispersed throughout the entire region. Less mobile species such as scallop and crab exhibited relatively more restricted areas of importance, which were associated with prevailing current direction, and proximity to shore in the case of green crab (Fig. 4). Broad areas of relatively high values for connectivity tend to occur when multiple sites are geographically proximate such as for cod populations on the westerns side of Newfoundland. Greater values for connectivity also occurred where land barriers constricted current pathways such as around Cape Breton Island, through the Strait of Canso, and through the Strait of Belle Isle.

Functional connectivity
The multiplication of the migration surfaces from EEMS by the current density maps of seascape connectivity effectively scales the importance of seascape connectivity values to gene flow. In particular, areas with lower-than-average migration from the EEMS rasters reduced the magnitude of seascape connectivity values within those areas, creating potential barriers not evident when considering connectivity through seascape alone (Fig. 5). This includes the gene flow barrier along the Scotian Shelf separating the northern and southern ecotypes and the isolation of the populations such as the Southern Newfoundland Scallop population. Like the conductance maps based on habitat resistances, areas including the nearshore waters off Cape Breton, and the Strait of Belle Isle between Newfoundland and Labrador, consistently showed areas of high functional connectivity, suggesting their potential importance in conservation planning.

Discussion
Connectivity has long been established as a core component of terrestrial conservation planning where areas of suitable (transmittable) habitat help to inform conservation design (Santiago and Pascual-Hortal 2007). To accomplish this, planners should ensure that conserved areas are sized, spaced, and positioned in a way that facilitates connectivity of organisms of conservation interest (Moffitt et al. 2011;Burgess et al. 2014). Here, we develop a novel method to highlight regions of high functional connectivity in the marine realm in a spatial context, directly integrating genetic structure with spatial transmittance as a product of both pelagic and benthic life phases. While previous studies in aquatic systems have investigated the correlation of genetic structure and resistance Fig. 4 Output from current density analysis using Circuitscape depicting effective resistance (inverse of density) for each of four species in this study pathways based on biophysical models (e.g., Thomas et al. 2015;Xuereb et al. 2018), few have directly integrated gene flow with dispersal potential into measures of spatial connectivity (e.g., Kininmonth et al. 2010). Measuring and implementing connectivity in marine conservation planning is challenged by the three-dimensional, transitory, aqueous environment in which animals move (Selkoe et al. 2016). For meroplanktonic species, processes regulating connectivity in pelagic early life history phases are distinct and often disconnected from the benthic processes that regulate connectivity later in life. Additionally, barriers to gene flow (e.g., mountains or landmasses) are often less conspicuous in the ocean and can reflect a combination of depth, currents, substrate, and localized environmental conditions ). This complexity challenges conventional conservation planning where conservation features and design elements are weighted spatially and incentivizes a different approach.
Integrating genetic and seascape connectivity into functional connectivity The number of approaches for incorporation of connectivity in marine conservation design has increased over the past two decades, for the most part motivated by high-level commitment for 'well-connected' networks of marine reserves (CBD 2010) and a strong theoretical basis relating resiliency of a network design with the relative degree of connectivity (e.g., Berumen et al. 2012). Simple approaches use minimum spacing to broadly categorize inferred connectivity via distance thresholds, but do not account for variation among taxa or within the planning region (spatial variation) (e.g., Friesen et al. 2019). More complex models apply particle tracking simulations to generalize intraspecific variation in dispersal range among species or dispersal phenotypes (e.g., Riginos et al. 2019). Others evaluate connectivity from a spatial habitat perspective where much like terrestrial approaches, connectivity is evaluated as the transmissibility of a species through a conditional habitat (Dickson et al. 2018). Similarly, biophysical and habitat models have been combined to evaluate successful dispersal and productive capacity (e.g., Magris et al. 2016), where dispersal simulations can explain spatial variation in gene flow (Pujolar et al. 2013;Xuereb et al. 2018). While each approach provides important information on connectivity at various scales relevant to conservation planning, there remain few approaches that integrate the biphasic dispersive Fig. 5 Maps of functional connectivity for all four species in this study, which are a product of the maps of gene flow (genetic connectivity) and resistance maps (structural connectivity) nature of many marine species and information on realized connectivity (gene flow) to represent connectivity spatially as is often applied in terrestrial settings (e.g., Stralberg et al. 2020). The integration of population genetic structure is a key consideration when attempting to characterize connectivity in any environment. All four species investigated here show fine-scale population structuring and significant isolation by distance (IBD) and isolation by resistance (IBR) relationships, indicating that their populations are more genetically divergent the further they are apart (both geographically and based on realized constraints to dispersal). Atlantic cod showed the strongest IBR relationship of F ST with its adult habitat suitability model, while both Atlantic scallop and green crab were more closely associated with composite oceanography and adult habitat suitability layers, suggesting potentially more complex IBR relationships. Lobster showed the strongest correlation of F ST with great circle distance (i.e., a simple IBD relationship), though this species exhibits the overall weakest population structuring (Benestan et al. 2015) among species tested. Reinforced by a recently published study using copy number variants, Dorant et al. (2020) detected very weak population structure among 21 populations of American lobster, reflecting high gene flow and large effective population sizes in this species. In each of the four species investigated here, habitat suitability shows a relatively continuous distribution across the latitudinal range with no conspicuous barriers to dispersal aside from landforms. Evaluating the distribution of potential connectivity through habitat using Circuitscape, in absence of any other information detailing potential barriers to functional connectivity, would likely overestimate connectivity between sites on either side of the cryptic genetic break described in Stanley et al. (2018). For example, in Atlantic Scallop, it was demonstrated that scaling habitat models based on this genetic break improved model performance, suggesting that the spatial structure previously identified was driven by environmental adaptation and would provide important context to evaluating habitat-based connectivity across the study area (Lowen et al. 2019). The use of EEMS provides a valuable estimate of gene flow through the planning area and thus accounts for reproductive success of migrants; however, it assumes equal ease of transmission throughout the seascape and reflects long-term biological processes (genetic structure). The combination of genetic information like EEMS with spatial information reflecting shorter term dynamics, like habitat suitability, should be cautioned when an IBR relationship is not present, which suggests another, unaccounted for process is driving (or has driven) genetic structure. However, when IBR is present, the incorporation of estimated effective migration values with seascape connectivity metrics can account for realized patterns of gene flow, intensifying connectivity values in areas with high estimated gene flow while downscaling connectivity values in areas with low gene flow, often revealing barriers to connectivity (e.g., Bertola et al. 2020).
As a basis for our approach, we applied previously developed species distribution models to evaluate potential pathways for connectivity across the study area. However, these models reflect the distribution of each species' post-pelagic life history and thus do not account for the pelagic dispersal component of these species' life histories. To address this in our framework, we also developed a method to assess potential pathways for larval dispersal using prevailing currents in the region. This approach was particularly important for species with low dispersal rates during their benthic life history. For example, both Browns Bank and Georges Bank are suitable habitat for Atlantic Scallop, while the Fundian Channel dividing these banks is not (Lowen et al. 2019). Gene flow and recruitment between these productive areas is unlikely to be achieved solely by successive migration and recruitment around the Fundian Channel and through the Gulf of Maine (100 km), when their planktonic larvae can simply disperse the relatively short distance (~ 60 km) between banks through the water column. While individual based hydrodynamic models and higher temporal resolution circulation model outputs could be used to infer pelagic dispersal and determine potential connectivity among conservation areas (e.g., Thomas et al. 2015), we sought to develop a first order estimate of resistance using prevalent currents (sensu Dambach et al. 2016) and dispersal bearing that could produce a raster of resistance compatible with Circuitscape. These raster layers could then be used to directly compare and integrate with habitat-based resistance maps to characterize pathways of least resistance and potential planktonic or pelagic dispersal among sampling sites.
Regions that showed consistently high levels of functional connectivity among these four species include southwestern Nova Scotia, northern Cape Breton, and the Strait of Belle Isle. The waters around southwest Nova Scotia serve as an area of both suitable habitat and high gene flow among populations in the Gulf of Maine, Bay of Fundy, and the western Scotian Shelf. Northern Cape Breton is a region of high connectivity between the Gulf of St. Lawrence and eastern Scotian Shelf and is an important migratory passage for a variety of marine species, including sharks, cetaceans, and sea turtles (e.g., O'Boyle 2012). Similarly, the Strait of Belle Isle serves as a region of high functional connectivity between the Labrador Sea, northern Newfoundland, southern Newfoundland, and the Gulf of St. Lawrence and is a known passage used by species such as Atlantic Cod and Atlantic Salmon as they move into or out of the Gulf of St. Lawrence during seasonal migrations, or to access feeding grounds (Templeman 1979;Strøm et al. 2017). Movement patterns can provide information for confirming areas of high functional connectivity outlined in our analysis. For example, detections of acoustically tagged Atlantic cod transiting through the St. Anns Bank Marine Protected Area north of Cape Breton (see Stanley et al. 2016 for details on the acoustic array) corresponds to an area of high functional connectivity identified for cod in this analysis. Deployment of acoustically tagged animals and receivers across a gradient of connectivitybased resistances outlined in our approach could further work to test the outcomes of mapping functional connectivity across a planning area.

Future directions
While the goal of this study was to develop a reproducible method for revealing regions of functional connectivity in the marine realm for use in marine spatial and conservation planning, populations were not originally sampled with this explicit intent. All samples for analyses of population genetics were collected in previous studies and were not collected in Marine Protected Areas or with marine conservation planning in mind. We emphasize that we have developed a method here to spatially visualize functional connectivity in the ocean, and that in order to apply these data layers targeted sampling from proposed conservation areas, spawning banks, nurseries, or known stocks will be needed, depending on the life stages and species desired to be protected. Ultimately, the sampling strategy will impact both the genetic interpolations and extent of functional connectivity, and so representative or extensive sampling sites across the geographic area of interest are needed (though EEMS is theoretically robust to sampling strategy if populations are sampled on both sides of the gene flow barriers; Petkova et al. (2016)). Further, as with most research studying landscape or seascape genetics, we used an aggregated population-based sampling scheme, whereas an individual-based sampling strategy may be more appropriate to sample over a broader spatial scale that corresponds with finer scale habitat heterogeneity (Prunier et al. 2013). Depending on the goal of a particular study, aggregate (e.g., F ST or other metrics) or individual-based sampling (e.g., in STRU CTU RE, Pritchard et al. 2000) may show advantages or disadvantages, though both methods should be applicable provided the genetic results can be interpolated across space and the habitat suitability models overlap with all sampled populations or individuals. In the coastal environment for example, impermeable landmasses can result in constraints on movement, resulting in constrained dispersal through passages such as the Strait of Canso and the Strait of Belle Isle. For species that are depth-limited and thus restricted to shore, including the intertidal zone, the use of habitat models alone can result in barriers when conducted at large spatial scales, which can be problematic when investigating connectivity of species with no pelagic larval stages. It is important that the resolution of genetic sampling, habitat suitability, and oceanographic models also need to reflect the study's objectives.
Though we have demonstrated the utility and integration of EEMs and Circuitscape outputs, there are a variety of other software tools that could potentially be integrated into this workflow. These include Omniscape (Landau et al. 2021), which implements circuit theory omnidirectionally across the seascape independent of user-defined nodes. Coalescent-based approaches like MAPS (migration and populationsize surfaces; Al-Asadi et al. 2019) provide a similar spatial representation of gene flow but with different underlying assumptions and data. Importantly our integrative framework can be adapted as new analytical approaches become available.

Applications in conservation network design
In Canada, a marine conservation target of 10% protected marine and coastal areas was achieved in 2020 and a new target of 30% by 2030 is now well underway. These marine conservation areas have the overall objective of conserving aquatic biodiversity through representative and well-connected protected areas (Schram et al. 2019). While individual sites have been developed or announced, and a national network is currently being developed, metrics of connectivity were not explicitly considered in regional network designs (Balbar et al. 2020). The results of our method produce spatial layers of functional connectivity that reveal important regions for connectivity, which can be imported directly into spatial conservation planning software such as Marxan (e.g., by setting an explicit spatial target for habitat weighted by its value for connectivity for a specific species or aggregate of species, by summing across functional connectivity layers). For example, our analysis can prioritize intermediate habitat between areas of focus (e.g., between known aggregations, populations, or existing protected areas) relevant for conservation priority species within a conservation network. Similarly, areas of high overlap between species (as shown for the Strait of Belle Isle in this study) can help focus conservation efforts towards those areas that support functional connectivity within a planning region. By representing connectivity spatially, these layers can be assigned targets and incorporated into traditional systematic conservation planning exercises. By overlaying existing networks of protected areas on species or aggregate functional connectivity maps, managers can identify gaps (gap analysis) within a network design that could be addressed with the addition of new sites (e.g., 'stepping-stone' locations) or enhanced conservation measures outside of MPAs. In this way, our approach provides spatial connectivity information to identify locations that serve as dispersal links through unprotected habitats, and thus can be used to augment existing (or draft) networks of MPAs to better promote connectivity overall. This method can also be revisited at regular intervals to reassess conservation area networks design over time or monitor conservation priorities as genetically divergent populations shift their range in response to climate change. For example, Lehnert et al. (2018) calculated a significant southward shift in the centre of the geographic cline of green crab ecotypes over a period of 15 years, driven by ongoing hybridization between northern and southern ecotypes. It may be prudent to conduct regular updates on a decadal scale for connectivity maps developed as species ranges shift and the global climate changes.
In summary, we have developed a method that highlights areas of seascape that are important for connectivity based on resistance derived from oceanography, adult habitat, and genetics using circuit theory for marine organisms to aid in conservation planning. The workflow developed here is broadly applicable in marine settings where genetic, oceanographic, and distribution models are available, and can be refined over time with new or higher-resolution genetic and spatial data. The resulting spatial connectivity information provides another feature to be considered among a suite of other criteria to inform the establishment and validation of MPAs and MPA networks.

Supporting information
Example output from EEMS (Online Resource 1), seasonality of pelagic larvae (Online Resource 2), correlations between linearized F ST and effective resistance distances (Online Resource 3), and F ST tables by species (Online Resources 4-7) are available online. The authors are solely responsible for the content and functionality of these materials. Queries (other than absence of the material) should be directly sent to the corresponding author.
Healthy Oceans Network and its partners: Department of Fisheries and Oceans Canada and INREST (representing the Port of Sept-Îles and City of Sept-Îles).
Data availability Not applicable.

Declarations
Competing interests The authors have no relevant financial or non-financial interests to disclose.

Ethical approval Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.