Seascape Genomics Reveals Metapopulation Connectivity Network of Paramuricea biscaya in the Northern Gulf of Mexico

The degree of connectivity among populations influences their ability to respond to natural and anthropogenic stressors. In marine systems, determining the scale, rate, and directionality of larval dispersal is therefore, central to understanding how coral metapopulations are interconnected and the degree of resiliency in the event of a localized disturbance. Understanding these source-sink dynamics is essential to guide restoration efforts and for the study of ecology and evolution in the ocean. The patterns and mechanisms of connectivity in the deep-sea (>200 m deep) are largely understudied. In this study, we investigated the spatial diversity patterns and metapopulation connectivity of the octocoral Paramuricea biscaya throughout the northern Gulf of Mexico (GoM). Paramuricea biscaya is one of the most abundant corals on the lower continental slope (between 1,200 and 2,500 m) in the GoM. The 2010 Deepwater Horizon oil spill (DWH) directly impacted populations of this species and thus are considered primary targets for restoration. We used a combination of seascape genomic analyses, high-resolution ocean circulation modeling, and larval dispersal simulations to quantify the degree of population structuring and connectivity among P. biscaya populations. Evidence supports the hypotheses that the genetic diversity of P. biscaya is structured by depth, and that larval dispersal among connected populations is asymmetric due to dominant ocean circulation patterns. Our results suggest that there are intermediate unsampled populations in the central GoM that serve as stepping stones for dispersal. The data suggest that the DeSoto Canyon area, and possibly the West Florida Escarpment, critically act as sources of larvae for areas impacted by the DWH oil spill in the Mississippi Canyon. This work illustrates that the management of deep-sea marine protected areas should incorporate knowledge of connectivity networks and depth-dependent processes throughout the water column.

The degree of connectivity among populations influences their ability to respond to natural and anthropogenic stressors. In marine systems, determining the scale, rate, and directionality of larval dispersal is therefore, central to understanding how coral metapopulations are interconnected and the degree of resiliency in the event of a localized disturbance. Understanding these source-sink dynamics is essential to guide restoration efforts and for the study of ecology and evolution in the ocean. The patterns and mechanisms of connectivity in the deep-sea (>200 m deep) are largely understudied. In this study, we investigated the spatial diversity patterns and metapopulation connectivity of the octocoral Paramuricea biscaya throughout the northern Gulf of Mexico (GoM). Paramuricea biscaya is one of the most abundant corals on the lower continental slope (between 1,200 and 2,500 m) in the GoM. The 2010 Deepwater Horizon oil spill (DWH) directly impacted populations of this species and thus are considered primary targets for restoration. We used a combination of seascape genomic analyses, high-resolution ocean circulation modeling, and larval dispersal simulations to quantify the degree of population structuring and connectivity among P. biscaya populations. Evidence supports the hypotheses that the genetic diversity of P. biscaya is structured by depth, and that larval dispersal among connected populations is asymmetric due to dominant ocean circulation patterns. Our results suggest that there are intermediate unsampled populations in the central GoM that serve as stepping stones for dispersal. The data suggest that the DeSoto Canyon area, and possibly the West Florida Escarpment, critically act as sources of larvae for areas impacted by the DWH oil spill in the Mississippi Canyon. This work illustrates that the management of deep-sea marine protected areas should incorporate knowledge of connectivity networks and depth-dependent processes throughout the water column.
Keywords: population genomics, connectivity, Gulf of Mexico, coral, seascape genomics, RADseq, dispersal, deep-sea INTRODUCTION Marine ecosystems have traditionally been considered "open" with few apparent barriers to dispersal. However, phylogeographic studies often reveal unexpected levels of population structuring or even previously unrecognized cases of cryptic speciation (Hellberg, 2009;Hoffman et al., 2012;Cerca et al., 2021). These studies have primarily focused on coastal ecosystems and species of significant economic importance. In comparison, the patterns and mechanisms that generate genetic diversity in the deep-sea (> 200 m deep) are largely understudied (Baco et al., 2016;Taylor and Roterman, 2017).
One general pattern in the deep-sea is that populations found at different depths (vertically separated by tens to hundreds of meters) are generally more differentiated than populations found at similar depths over large geographical areas (horizontally separated by hundreds to thousands of kilometers) (Taylor and Roterman, 2017). However, the mechanisms responsible for this pattern remain poorly understood. Determining the scales of connectivity of marine populations and the mechanisms behind them is crucial for the conservation of marine ecosystems (Palumbi, 2003;Kinlan et al., 2005;Botsford et al., 2009;Gaines et al., 2010), and the study of diversification and evolution in the ocean (McClain and Mincks Hardy, 2010).
Population genetic methods enable the identification of genetic structuring patterns and estimate the scale, rate, and direction of reproductive exchange among marine populations (Breusing et al., 2016;Galaska et al., 2017;Bertola et al., 2020). These inferences, when coupled with analyses of environmental parameters, physical models of ocean circulation, and simulations of larval dispersal, can significantly enhance our understanding of connectivity networks at scales relevant to management (Benestan et al., 2016;Sandoval-Castillo et al., 2018;Xuereb et al., 2018;Bernatchez et al., 2019;Bracco et al., 2019). This integrative approach is known as seascape genetics (Galindo et al., 2006;Selkoe et al., 2016). Only a handful of studies have implemented seascape approaches in the deep sea. These have predicted the presence of intermediate "phantom" populations of hydrothermal vent species along mid-ocean ridges (Breusing et al., 2016) and have suggested that variables related to currents and food sources may explain a significant fraction of observed genetic patterns of sponge and coral species (Zeng et al., 2020).
Corals are essential foundational species in deep-sea benthic habitats and are typically slow-growing and long-lived (Roark et al., 2009;Sherwood and Edinger, 2009;Prouty et al., 2011Prouty et al., , 2016Girard et al., 2019). Deep-sea coral ecosystems are analogous to islands in that they are discrete and spatially separated. Each community serves as an oasis or biodiversity hotspot by locally enhancing the abundance and diversity of invertebrates and fishes (Henry and Roberts, 2007;Ross and Quattrini, 2007;Cordes et al., 2008;Rowden et al., 2010;Demopoulos et al., 2014). These characteristics of deep-sea corals make them particularly susceptible to anthropogenic impacts and a priority for conservation efforts.
The degree of connectivity among deep-sea coral populations influences the probability of speciation (Quattrini et al., 2015;Herrera and Shank, 2016) and likely contributes to their ability to respond to natural and anthropogenic stressors. Determining the scale, rate, and directionality of larval dispersal is therefore, central to understanding how coral metapopulations are interconnected and the degree of resiliency in the event of a localized disturbance, such as an oil spill (Jones et al., 2007;Almany et al., 2009). Understanding these source-sink dynamics is essential to guide restoration efforts (Lipcius et al., 2008;Puckett and Eggleston, 2016).
Herein, we investigate the spatial patterns of genetic variation and metapopulation connectivity of the octocoral Paramuricea biscaya throughout the northern Gulf of Mexico (GoM), using a seascape genomics framework. Paramuricea biscaya, distributed between 1,200 and 2,500 m, is one of the most common and abundant corals on hardgrounds on the lower continental slopein the GoM . The 2010 Deepwater Horizon oil spill (DWH) directly impacted populations of this species (White et al., 2012;Fisher et al., 2014;DeLeo et al., 2018) and thus are considered primary targets for restoration (Deepwater Horizon Natural Resource Damage Assessment Trustees, 2016). We use a combination of population seascape genomic analyses, high-resolution ocean circulation modeling, and larval dispersal simulations to quantify the degree of structuring and connectivity among DWH impacted and non-impacted populations. This paper is a companion to the paper by Liu et al. (2021) that describes the ocean circulation modeling and larval dispersal simulations. Here we test the hypothesis that the genetic diversity of P. biscaya is predominantly structured by depth, and to a lesser degree, by distance. We also test the hypothesis that larval dispersal among connected populations is asymmetric due to dominant ocean circulation patterns.

Collection of Samples
We sampled Paramuricea biscaya colonies from six sites in the Northern Gulf of Mexico at depths between 1,371 and 2,400 m ( Table 1 and Figure 1). The 2010 Deepwater Horizon oil spill directly impacted P. biscaya populations at three of these sites in the Mississippi Canyon area (MC294, MC297, and MC344) (White et al., 2012;Fisher et al., 2014). . We imaged individual coral colonies before and after removing a small distal branch using hydraulic manipulations mounted on remotely operated vehicles or submarines. We stored samples in insulated containers until the recovery of the vehicles by the surface vessel. Subsamples of each specimen were preserved in liquid nitrogen or 95% ethanol and stored at −80 • C.

Molecular Laboratory Methods
To characterize the genetic diversity of P. biscaya individuals, we performed reduced representation DNA sequencing (RADseq) (Baird et al., 2008;Reitzel et al., 2013). DNA was  (Nei, 1978). See the main text for the source of the other environmental parameters. Environmental parameters represent averages (and standard deviations) of values extracted from gridded datasets using geographical coordinates of each sample.
purified using the Qiagen DNeasy Blood and Tissue Kit following manufacturers' protocols. We checked DNA integrity and purity by visual inspection on a 1% agarose gel and a NanoDrop spectrophotometer (NanoDrop Technologies), respectively. DNA concentration was determined and normalized using a Qubit 4.0 fluorometer (Invitrogen). We confirmed species identification through DNA barcoding of the COI mitochondrial gene following the protocols described by Quattrini et al. (2014) (NCBI GenBank Accession numbers MT795490 to MT795554). Floragenex Inc. (Eugene, OR) performed RAD sequencing library preparation utilizing the 6-cutter PstI restriction enzyme on quality-checked and concentration-normalized high-molecularweight DNA. Using the program PredRAD (Herrera et al., 2015), we predicted tens of thousands of cleavage sites in coral genomes with the PstI restriction enzyme. Libraries were dual-barcoded and sequenced on an Illumina Hi-Seq 4,000 1 × 100 platform.

Data QC and Single Nucleotide Polymorphism Calling
We de-multiplexed and quality filtered raw sequence RADseq reads using the process_radtags program in Stacks v2.1 (Catchen et al., 2013) with the following flags: -inline_null, -r, -c, and -q, with default values. We performed read clustering and single nucleotide polymorphism (SNP) calling using the DeNovoGBS (Parra-Salazar et al., 2021) module of the software package NGSEP v4.0.1 (Tello et al., 2019). This software is more computationally efficient and has comparable or better accuracy than programs like Stacks or pyRAD (Eaton, 2014) for de novo analysis of genotype-by-sequencing data (Parra-Salazar et al., 2021). We assumed a heterozygosity rate of 1.5% (-h 0.015) as calculated from the short read genome-wide data of the sister species Paramuricea sp. type B3  using the software GenomeScope v2.0 (Vurture et al., 2017) from the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under BioProject number PRJNA574146 (Vohsen et al., 2020).
The all_snp dataset was imported into the R v4.0.3 statistical environment (R Core Team, 2013) for further filtering. We excluded individuals if they had missing data in more than 35% of the SNP loci or identified as clones by clonecorrect function from the R package poppr v2.8.6 (Kamvar et al., 2014). We excluded SNP loci if their observed heterozygosity was greater than 0.5, as estimated with hierfstat v0.5 , or if their allelic frequencies were not in Hardy-Weinberg equilibrium, as estimated with pegas v0.14 (Paradis, 2010) (B = 1,000, p < 0.01). We randomly retained one SNP per RAD locus to reduce the risk of violating the assumption of independence among SNP. Finally, 10 SNPs in RAD loci identified as potentially under positive selection by BayeScan were excluded (Supplementary Figure 1). This dataset, containing 4,248 unlinked neutral SNPs across 133 individuals, is hereafter referred to as the neutral dataset.

Genetic Connectivity
To measure the genetic connectivity among sampling sites, we estimated migration rates (m), defined as the proportion of immigrant individuals in the last two generations, using BAYESASS v3.0.4.2 (Wilson and Rannala, 2003). Twelve independent runs with different random seeds were performed using the neutral dataset. We ran each analysis for 100 million Markov chain Monte Carlo (MCMC) iterations, with 50 million burn-in iterations and one thousand iterations sampling frequency. Mixing parameters (-m0.35 -a0.9 -f0.09) were optimized to ensure adequate mixing (acceptance rates between 20 and 60%). MCMC trace files were examined in the program Tracer v1.7.1 (Rambaut et al., 2018) to evaluate convergence and consistency of estimates among runs. We calculated point estimates of m as the median of the posterior Frontiers in Marine Science | www.frontiersin.org distribution and their uncertainty as 95% High Posterior Density (HPD) intervals.

Potential Connectivity
To identify dispersal mechanisms that could explain genetic connectivity estimates, we compared our results with the potential connectivity estimates (probability of connectivity through larval dispersal among sampling sites) by Liu et al. (2021). Briefly, Liu et al. (2021) simulated the dispersal trajectories of neutrally-buoyant Lagrangian particles in an implementation of a high-resolution three-dimensional Coastal and Regional Ocean COmmunity hydrodynamic model (CROCO) (Liu et al., 2021). The model encompassed the area between 98-82 • W and 24-31 • N and had a horizontal grid resolution of approximately 1 km and 50 vertical sigma (density) layers. Once per season, Lagrangian particles were deployed uniformly at the seafloor in 0.05 • × 0.05 • boxes (approximately 5 km by 5 km) centered at the location of the sampling sites, with a zonal (dx) and meridional (dy) intervals of 0.00075 • (approximately 75 m). This configuration resulted in 4,489 particles released at each site, and nearly 27,000 particles in each seasonal release. The particles were tracked offline using the Lagrangian tool Ichthyop (Lett et al., 2008) and recorded hourly. Horizontal connectivity through larval dispersal among sampling sites (l h ) was defined as the average proportion of neutrality-buoyant Lagrangian particles released at a source site (i) area that passed over another site (j) area (sink) after 56 days (computational constraints limited the length of the tracking) starting from January 25th, April 25th, July 24th, and November 1st, 2015. The pelagic larval duration (PLD) for Paramuricea biscaya is unknown, but Hilario et al. (2015) found that a PLD between 35 and 69 days seems representative of 50-75% of deep-sea species. The definition of vertical connectivity (l v ) is the same as horizontal connectivity, except that a particle also has to pass within 50 m of the sink site's seafloor depth. Liu et al. (2021) also evaluated longer PLDs by extending the Lagrangian tracking starting November 1st to 148 days and with additional Eulerian dye releases followed for 120 days. The dye release indicated that although Lagrangian particles cover a smaller area than the Eulerian dye, they capture the same main dispersal features and do not predict substantially different connectivity patterns. No other biological parameters such as larval growth, mortality, settlement, and swimming because they are unknown for the study species.

Population Genetic Structure
To determine the patterns of genetic structuring of the sampled P. biscaya corals, we performed a discriminant analysis of principal components (DAPC) on the neutral dataset using the R package adegenet v2.1.3 (Jombart, 2008). DAPC was performed with and without sampling locations as priors after estimating the optimal number of principal components with the function optim.a.score. For the DAPC with no priors, we applied the Bayesian Information Criterion to choose the optimal number of clusters (K) that explain the genetic variability in the dataset using the function find.clusters.
We also inferred population structuring patterns (as historical lineages) with the neutral dataset by maximizing the posterior probability of the genotypic data, given a set number of clusters (K). This method is known as Bayesian population clustering and is implemented in the program Structure v2.3.4 (Pritchard et al., 2000). We used the admixture model with uncorrelated allele frequencies. The MCMC was run for 1.1 × 10 6 repetitions (burnin period 1 × 10 5 ). We evaluated values for K from 1 to 6 (10 replicates each). We selected the optimal value of K using the program StructureHarvester v0.6.92 (Earl and vonHoldt, 2012) according to the ad hoc K statistic (Evanno et al., 2005), which is the second-order rate of change of the likelihood function. We visualized structure results using the program R package starmie (Tonkin-Hill and Lee, 2016).
We performed a hierarchical Analysis of Molecular Variance (AMOVA) (Excoffier et al., 1992) with the neutral dataset to calculate F-statistics and test for differentiation at the individual, site, and genetic cluster levels. The AMOVA, performed in genodive v3.04 (Meirmans, 2020), assumed an infinite-alleles model. We calculated pairwise F ST (Weir and Cockerham, 1984) differentiation statistics among sampling sites with the R package assigner v.0.5.8 (Gosselin et al., 2020).

Redundancy Analyses
To quantify environmental variables' significance and relative importance in shaping genetic diversity in P. biscaya, we used a series of redundancy analyses (RDA) in the R package vegan v2.5 (Oksanen et al., 2007). RDA has two steps. First, a multiple linear regression between genetic (response) and environmental (explanatory) data matrices produces a matrix of fitted values. Second, a principal components analysis (PCA) of the fitted values. The PCA axes are linear combinations of the explanatory variables (Legendre and Legendre, 2012).
We performed site-level RDA (Legendre and Legendre, 2012) on sites' allelic frequencies and geographical distances. We first transformed geographical distances as in-water distances using the lc.dist function of the R package marmap v1.0.4 (Pante and Simon-Bouhet, 2013), and later represented as distance-based Moran's eigenvector maps (dbMEM) (Dray et al., 2006) using the R package adespatial (Dray et al., 2018).
We bottom seawater potential density (σ θ ) values using the R package oce v1.2 (Kelley, 2018). We obtained average monthly surface chlorophyll concentration and primary productivity values between 2011 and 2018 from the E.U. Copernicus Marine Service Information, Copernicus Globcolour ocean products grids OCEANCOLOUR_GLO_CHL_L4_REP_OBSERVATIONS _009_082 and OCEANCOLOUR_GLO_OPTICS_L4_REP_OBS ERVATIxONS_009_081 (accessed on June 2021), 1 and summarized as mean and standard deviation grids with a 4 km resolution. Individual parameter values were extracted from these grids using the latitude and longitude of each sampled coral.
To avoid problems with highly correlated environmental variables (Dormann et al., 2013), we performed a pairwise correlation test and removed variables with a correlation coefficient |r| > 0.7 and a p-value < 0.05. We retained the most seemingly ecologically relevant variable when two or more variables were correlated. We evaluated the explanatory importance of each environmental variable using forward selection and analysis of variance (ANOVA) after 10,000 permutations (α = 0.05) using the ordistep function in vegan. Retained environmental variables were included in the dbRDA using the dbrda function in vegan. We performed a variance partitioning analysis with the function varpart and tested its significance through global and marginal ANOVAs (1,000 permutations, α = 0.01).

Population Connectivity
Migration rates (m) among sites estimated from genetic data using BAYESASS were overall low (average = 0.011, standard deviation = 0.065), with a few exceptions. Approximately

Population Genetic Structure
DAPC analysis with no location priors indicated that there is metapopulation substructuring within P. biscaya's sampled range. The variability in the genetic data was explained by two clusters K D1 (n = 89) and K D2 (n = 44) (optimal K = 2, BIC = 616.8, 18 retained PCs; Figure 1C). Each individual was considered a member of the group with the highest probability. The first cluster is mainly composed of individuals collected at sites DC673 (100% of sampled individuals belong to K D1 ), KC405 (100% K D1 ), and MC344 (83% K D1 ), while the second cluster is mainly composed of individuals collected at sites MC294 (93% of sampled individuals belong to K D2 ), MC297 (54% K D2 ), and GC852 (75% K D2 ) (Figures 1C-E). The first discriminant axis, calculated by DAPC analysis with location priors, explained 58.4% of the variance and primarily reflected the differentiation between the two inferred clusters K D1 and K D2 . This differentiation seemed to be associated with depth as samples assigned to K D1 were on average found at deeper locations (mean depth at which K D1 individuals were sampled: x KD1 = 1,896 m, standard deviation: s KD1 = 308 m) than individuals assigned to K D2 (x KD2 = 1454 m, s KD2 = 130 m) ( Figure 1F).
The STRUCTURE analyses of Bayesian population clustering confirmed the presence of two ancestry clusters, K S1 (n = 91) and K S2 (n = 42), that largely corresponded to the clusters identified by the DAPC K D1 and K D2 , respectively (we considered each individual a member of the group for which it had the highest membership probability Q). To maintain consistency with other studies in corals (Carlon and Lippé, 2011;Serrano et al., 2016), we defined an admixed individual as having a Q > 0.1 for both clusters. These analyses indicate that, overall, 15% of individuals have an admixed ancestry (Figure 1G), but proportionally there are more admixed individuals assigned to K S2 (24% of individuals) than to K S1 (11%). Ancestry cluster K S1 is dominant in DC673 (mean probability of membership at that site:x QS1 = 0.97), KC405 (x QS1 = 0.98) and MC344 (x QS1 = 0.79), whereas K S2 is dominant in MC294 (x QS2 = 0.88), GC852 (x QS2 = 0.65), and MC297 (x QS2 = 0.52). Pairwise F ST (Weir and Cockerham, 1984) statistics among sites were consistent with the DAPC and STRUCTURE results showing greatest differentiation among sites with a majority of individuals assigned to different K clusters, and lowest among sites with a majority of individuals assigned to the same K cluster ( Figure 1H). The AMOVA analysis indicated that 11.4% of the observed genetic variation could be attributed to differences among individuals (F IS = 0.117, p = 0.001), 0.3% to differences among sites (F SC = 0.003, p = 0.001) and 1.8% to differences between DAPC clusters K D1 and K D2 (F CT = 0.018, p = 0.001).

Redundancy Analyses
Site-level RDA failed to detect a significant correlation (α = 0.05) between geographic distance (as dbMEM eigenvectors) and genetic differentiation, thus rejecting the hypothesis of isolation by distance.
Of the environmental variables, we excluded bottom temperature and salinity from the individual-level dbRDA as the ranges of these parameters across the study sites were too small to be biologically important (4.21-4.84 • C, and 34.95 and 35.10 PSU). Mean surface primary productivity (retained for dbRDA) was significantly correlated with its standard deviation and latitude and mean and standard deviation of surface chlorophyll concentration. Depth (retained) was significantly correlated with bottom seawater potential density (σ θ ). Mean bottom oxygen concentration (retained), longitude (retained), and mean bottom current speed (retained) were not significantly correlated with any other environmental variable.
We incorporated depth, mean bottom dissolved oxygen concentration, mean surface primary productivity, longitude, and mean bottom current speed into an initial dbRDA model as these were the only significant independent variables identified by forward selection (ANOVA, p-values < 0.05). These variables significantly contributed to the model, except for mean bottom current speed (ANOVA, p-value > 0.01), which was subsequently excluded.
Globally, the percentage of the genetic variation explained by environmental variables was 7.37% ( Table 2). Depth [collinear with bottom seawater potential density (σ θ )] had the largest effect (explaining 3.8% of the variance), followed by mean bottom oxygen concentration (2.0%), mean surface primary productivity (collinear with six other variables, see above) (1.8%), and longitude (1.1%). The combined effect of depth and density is evident in the dbRDA plots (Figure 3). dbRDA axes are linear

DISCUSSION
Population Connectivity: Scale, Rate, and Directionality Larval dispersal simulations in the study area show a prevailing westward pathway of dispersal along isobaths in the 1,000-2,000 m range in all seasons (Liu et al., 2021). Long-distance dispersal (more than 100 km) driven by strong deep recirculation currents (Bracco et al., 2016) may occur for larvae originating in the DeSoto Canyon area (DC673) (Liu et al., 2021). These larvae can reach the Mississippi Canyon area in less than 2 months (Liu et al., 2021), thus explaining the source-sink dynamics identified between these sites by migration rate estimates (m) from genetic data (Figure 2). These source-sink dynamics are highly depthdependent. Our estimates suggest that 15-27% of individuals at MC344 (1,852 m) likely immigrated from the De Soto Canyon area (DC673, 2,254 m) within the last one or two generations. For MC297 (1,577 m), 3-15% are likely immigrants from DC673, and less than 8% for MC294 (1,371 m). The limiting effect of depth on vertical connectivity is most striking within the Mississippi Canyon (Figure 2), where we found no evidence of substantial gene flow among sites. The limited amount of vertical diapycnal mixing possible over the short horizontal distances that separate them (tens of kilometers, see the following section for further discussion on the role of depth) may explain the limited 3D connectivity among these sites (Bracco et al., 2019;Liu et al., 2021).
Our analyses indicate that the population of P. biscaya at DC673 should be a conservation priority to restore the impacted populations at MC344 and MC297. Additional sampling and modeling throughout P. biscaya's depth range (1,000-2,600 m) in the DeSoto Canyon and West Florida Escarpment are necessary to fully understand the role of this region as a source of larvae for DWH impacted populations in the Mississippi Canyon and identify other sites in need of protection.
Larval dispersal models predict that larvae originating from the Keathley Canyon area (KC405) can disperse the furthest (maximum horizontal distances 154 and 426 km after 56 and 148 days, respectively, Liu et al., 2021). The highly variable currents that characterize this area can explain this potential for long-distance dispersal (Liu et al., 2021). However, these models fail to predict the degree of direct genetic connectivity estimated between KC405 and DC673 (Figure 2). Similarly, the relative importance of the Green Canyon site GC852 as a source of larvae to the Mississippi Canyon area, indicated by the migration rate estimates, is not consistent with the connectivity probabilities estimated by the numerical larval simulations (Liu et al., 2021). The dispersal distances for larvae out of GC852 do not seem to exceed 100 km after 56 days (400 km after 148 days) (Liu et al., 2021). Thus no direct connectivity is predicted between the Green Canyon and Mississippi canyon sites separated by more than 300 km. Additional factors that may contribute to some of the inconsistencies between migration rates and connectivity probabilities include potential violations to the underlying assumptions of the larval dispersal (Liu et al., 2021), and the population genetic models. For example, the migration rate estimates calculated by BayesAss are most accurate and precise when true migration rates among populations are moderate (m < 0.333), population structuring is significant and sampling of individuals and loci is substantial (Faubet et al., 2007;Meirmans, 2014).
The patterns of genetic connectivity between KC405 and DC673, and GC852 and the Mississippi Canyon sites cannot be explained by larval dispersal modfels unless intermediate populations that act as stepping stones are included in the simulations (Liu et al., 2021) when the role of interannual variability is accounted for by using the advection pathways predicted by HYCOM data. Additional targeted exploration and sampling, informed by habitat suitability (Georgian et al., 2020) and dispersal models (Liu et al., 2021), are necessary to test this connectivity hypothesis and clarify the role of western populations in the restoration of DWH impacted populations.

Metapopulation Structuring by Depth
All of our analyses support the existence of two clusters or "stocks" of Paramuricea biscaya in our samples, both of which were impacted by the DWH oil spill. Previous studies sequenced mitochondrial DNA of P. biscaya (mtCOI + igr + MutS) and recovered three haplotypes of P. biscaya (B1, B1a, B2) in the northern Gulf of Mexico Quattrini et al., 2014;Radice et al., 2016). We found that mitochondrial haplotypes bear no direct correspondence with the genomic clusters (Supplementary Figure 3). Mitochondrial markers are well known for lacking sufficient variability at low taxonomic levels in octocorals and are subject to incomplete lineage sorting (Pante et al., 2015;Herrera and Shank, 2016;Quattrini et al., 2019). We suggest that mitochondrial DNA barcoding data should not be used to resolve differences at the population level, especially in the context of management and restoration, and in many cases at the species level, as it could lead to incorrect interpretations and inadequate policy decisions.
Geographic distance is not a significant variable structuring the genetic diversity of P. biscaya within the GoM. Despite only being separated by tens of kilometers, the populations in the Mississippi Canyon impacted by the DWH oil spill (MC294, MC297, and MC344) have distinct genetic compositions (Figure 1). The population's genetic composition at MC344 is most similar to those found at the DeSoto Canyon (DC673) and Keathley Canyon (KC405), hundreds of kilometers away. Consistent with results from previous studies of deep-sea populations (Taylor and Roterman, 2017), depth is a critical variable structuring the genetic diversity of P. biscaya. MC344 is the deepest of the three sites at the Mississippi Canyon (MC294: 1,371 m, MC297: 1,577 m, and MC344: 1,852 m), and its population is mainly composed of individuals whose ancestry is predominantly from the first cluster K D1 [83%; the other deep sites DC673 (2,254 m) and KC405 (1,679 m) are also almost entirely made up of individuals with K D1 ancestry]. MC294, the shallowest, has a population whose ancestry is mainly from the second cluster, K D2 (93%). MC297 sits at an intermediate depth and has a population of mixed ancestry, split roughly in half.
Seascape genomic analyses provide statistical support for the role of depth. The combined effect of depth and bottom seawater potential density (σ θ ) contributes the most toward explaining the genetic variability in P. biscaya among the environmental variables explored in this study (Figure 3 and Table 2). Among the environmental variables known to show collinearity with depth, hydrostatic pressure may be the most biologically important for Paramuricea biscaya. Hydrostatic pressure increases linearly with depth (at a rate of roughly 1 atmosphere every 10 m). Other variables, such as dissolved oxygen concentration, pH, temperature, and salinity, do not vary sufficiently within the depth and geographical range of the examined populations in the study area to exert any significant adaptive pressure that could drive diversification. Several studies have suggested that pressure can be a significant selective force in the deep sea, often driving the evolution of pressure-adapted enzymes and other biomolecules (Somero, 1992;Lan et al., 2017Lan et al., , 2018Gaither et al., 2018;Lemaire et al., 2018;Gan et al., 2020;Weber et al., 2020).
Another potentially important variable known to be collinear with depth is the flux of particulate organic matter from the surface ocean to the seafloor (POC flux). POC is the primary food source for most deep-sea organisms, and it is known to structure biodiversity patterns on the benthos (Woolley et al., 2016). POC flux decreases exponentially with depth Mouw et al., 2016) and could therefore be a significant selective force (Quattrini et al., 2017) and a major driver of biodiversity patterns in the deep sea. However, POC accumulates on the seafloor, where it can be resuspended through the interaction of bottom currents and complex topography (Wilson et al., 2015;Amaro et al., 2016). The role of POC resuspension is uncertain given the diversity of habitats where P. biscaya is found, from carbonate outcrops (e.g., GC852 and MC sites) to near-vertical walls (DC673, KC405). Furthermore, episodical delivery episodes of POC are challenging to incorporate in models but are likely biologically important (Smith et al., 2018). In situ measurements would be needed to quantify differences in food delivery at these sites. Although chlorophyll-a concentrations, sea surface temperature, and photosynthetically active radiation are used to model net primary productivity (Behrenfeld and Falkowski, 1997) and POC fluxes (Pace et al., 1987), their relationship is not always predictable. In the northern Gulf of Mexico, confounding factors such as planktonic community composition can cause discrepancies between modeled and in situ POC flux measurements (Biggs et al., 2008;Maiti et al., 2016).
Bottom seawater potential density (σ θ ) could play an important role if larvae behave as neutrally buoyant particles dispersing along isopycnals as suggested for other deep-sea corals and sponge species (Dullo et al., 2008;Kenchington et al., 2017;Bracco et al., 2019;Roberts et al., 2021). Larval dispersal along narrow density envelopes associated with water mass structuring may serve as a mechanism for increasing reproductive success and gene flow among deep-sea metapopulations of species with neutrally buoyant larvae while simultaneously facilitating prezygotic isolation by limiting dispersal across depth (Miller et al., 2011). The results from our potential connectivity analyses further suggest that the ocean circulation, and specifically the limited diapycnal mixing, may prevent neutrally-buoyant larvae from spreading across depth ranges.
Remarkably, all sampled sites, except for DC673, have a proportion of individuals with admixed ancestry, suggestive of successful crosses beyond F1 or F2 generations among clusters (Figure 1). This could be indicative of an absence of post-zygotic isolation barriers between clusters. There are two possible explanations for this pattern: incipient sympatric speciation or secondary contact. Incipient sympatric ecological speciation through niche specialization is a possible driver of the observed pattern of population structuring (González et al., 2018). Due to the relative environmental stability at the depth range of P. biscaya, it is plausible that specialization to pressure and food gradients would occur over the species' depth range and be reinforced by density-driven dispersal limitation. Alternatively, secondary contact could occur by recent colonization of the GoM by a P. biscaya lineage from the Caribbean Sea or the Atlantic Ocean. Differentiating between the two possibilities would require a combination of demographic modeling and additional sampling throughout the range of P. biscaya, i.e., not limited to the GoM.
The presence of these two genetic stocks should be taken into consideration for restoration activities that involve propagation in nurseries and transplantation (Baums et al., 2019). The possibility that the stocks are partially reproductively isolated and depth-adapted suggests that receiving populations would benefit most from transplants from populations to which they are already genetically connected.

CONCLUSION
In this study, we found support from population genomic analyses and larval dispersal modeling for the hypothesis that the P. biscaya metapopulation in the northern GoM is predominantly structured by depth, and to a lesser degree, by distance. Further, both lines of evidence (genetic and modeling) support the hypothesis that larval dispersal among connected populations is asymmetric due to dominant ocean circulation patterns. Utilizing a seascape genomic approach brought a more holistic understanding of the population connectivity of this species than either population genetics or modeling could on its own. There are likely intermediate unsampled populations that serve as stepping stones for dispersal. These may explain some of the observed genetic connectivity that could not be explained as direct dispersal by larval simulations. Although dispersal and connectivity patterns of organisms are highly species-dependent, the integrative framework of this study provides valuable insights to understand the connectivity of deep-sea metacommunities broadly (Mullineaux et al., 2018).
This study further illustrates that management of marine protected areas (MPAs) should incorporate connectivity networks and depth-dependent processes throughout the water column. Doing so could help preserve genetic diversity and increase species resilience to extreme climate events and anthropogenic impacts. We suggest that the DeSoto Canyon area, and possibly the West Florida Escarpment, critically act as sources of larvae that may repopulate areas impacted by the 2010 Deepwater Horizon oil spill in the Mississippi Canyon. Active management of these source sites is essential to the success of restoration efforts.