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). Collections took place during expeditions in 2009 (R/V Ron Brown, ROV Jason II), 2010 (R/V Ron Brown & R/V Atlantis, ROV Jason II & HOV Alvin), 2011 (MSV Holiday Chouest, ROV UHD-34), and 2017 (MSV Ocean Intervention II & MSV Ocean Project, ROV Global Explorer & ROV Comanche). 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.

To characterize the genetic diversity of P. biscaya individuals, we performed reduced representation DNA sequencing (RAD-seq) (Baird et al., 2008; Reitzel et al., 2013). DNA was 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-molecular-weight 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.

We de-multiplexed and quality filtered raw sequence RAD-seq 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 (Doughty et al., 2014) 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).

Single nucleotide polymorphisms were filtered using vcftools v0.1.16 (Danecek et al., 2011) to exclude SNP loci that: (1) had more than 30% missing data, (2) a mean depth of coverage smaller than 10x or greater than 100x, (3) a minor allele frequency smaller than 0.01, and (4) had more than two alleles. The resulting dataset, containing 12,948 SNPs and 154 individuals, is hereafter referred to as the all_snp dataset. BayeScan v.2.01 (Foll and Gaggiotti, 2008) was used to identify SNP potentially under positive selection (–n 5,000 –burn 50,000 –pr_odds 10,000, Qval < 0.05).

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 (Goudet et al., 2005), 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.

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 distribution and their uncertainty as 95% High Posterior Density (HPD) intervals.

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.

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⁶ repetitions (burn-in period 1 × 10⁵). 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).

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 performed individual-level distance-based RDA (dbRDA) (McArdle and Anderson, 2001) with the matrix of genetic distances calculated from the neutral dataset and a matrix of environmental variables. Missing genotypes in each individual were first imputed by assigning the most common genotype for each locus at the collection site. Environmental variables included: depth, latitude, longitude, bottom temperature, salinity, bottom current speed, bottom oxygen concentration, bottom seawater potential density (σ_θ), surface chlorophyll concentration, and surface primary productivity. Average monthly bottom temperature, salinity, current speed values between 2011 and 2018, as extracted by Goyert et al. (2021) from the HYbrid Coordinate Ocean Model (HYCOM) for the Gulf of Mexico, were summarized as mean and standard deviation grids with a 4 km resolution. Bottom oxygen concentration values are annual means gridded from the World Ocean Database by Goyert et al. (2021) at a 370 m resolution. We calculated 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),¹ 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).

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 20% of individuals at the sites MC344 (depth 1,852 m) and KC405 (1,679 m) likely immigrated from the De Soto Canyon area (site DC673, depth 2,254 m), within the last one or two generations [Figure 2A and Supplementary Figure 2; m_{DC673–MC344} = 0.210, 95% HPD_{DC673–MC344} = (0.145,0.269); m_{DC673–KC405} = 0.181, 95% HPD_{DC673–KC405} = (0.123,0.236)]. The potential contribution of immigrants from DC673 to site GC852 (depth 1,407 m), in the Green Canyon area, and sites MC294 (depth 1,407 m) and MC297 in the Mississippi Canyon area was smaller but still substantial, ranging between 3 and 10% [Figure 2A and Supplementary Figure 2; m_{DC673–GC852} = 0.097, 95% HPD_{DC673–GC852} = (0.050,0.149); m_{DC673–MC294} = 0.051, 95% HPD_{DC673–MC294} = (0.000,0.076); m_{DC673–MC297} = 0.085, 95% HPD_{DC673–MC297} = (0.028,0.153)]. These analyses also indicate that GC852 may also be an important source of immigrants to the Mississippi Canyon area. The potential contribution of immigrants from GC852 to sites MC344, MC297, and MC294 ranges between 5 and 24% [Figure 2 and Supplementary Figure 2; m_{GC852–MC344} = 0.052, 95% HPD_{DC673–GC852} = (0.012,0.105); m_{GC852–MC294} = 0.236, 95% HPD_{DC673–MC294} = (0.167,0.295); m_{DC673–MC297} = 0.139, 95% HPD_{DC673–MC297} = (0.071,0.212)].

Horizontal connectivity probabilities (l_h) calculated from larval dispersal simulations recovered a remarkable congruence with the estimated migration rates (m) concerning the role of the De Soto Canyon area DC673 as a source of larvae for the Mississippi Canyon sites (Figures 2B,D, l_{h DC673–MC344} = 0.184,; l_{h DC673–MCC294} = 0.051; l_{h DC673–MC297} = 0.096), but not for the Green Canyon (GC852) or the Keathley Canyon (KC405) areas. The congruence is only maintained in the 3D connectivity probability (l_v) for MC344 (Figures 2C,D, l_{v DC673–MC344} = 0.031). The larval dispersal simulations also predict bi-directional connectivity between MC294 and MC297 (l_{h MC294–MC297} = 0.348,; l_{h MC297–MCC294} = 0.172; l_{v MC294–MC297} = 0.082; l_{v MC297–MCC294} = 0.012) but the estimated migration rates between these sites are low [Figure 2A and Supplementary Figure 2; m_{MC294–MC297} = 0.013, 95% HPD_{MC294–MC297} = (0.000,0.051); m_{MC297–MCC294} = 0.012, 95% HPD_{MC297–MCC294} = (0.000,0.049)].

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).

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 combinations of the environmental variables. dbRDA axis 1, which explains 52.5% of the variation, broadly splits individuals belonging to different DAPC clusters. This differentiation is primarily driven by depth as indicated by the environmental variables vectors and suggested in Figures 1E,F.

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 depth-dependent. 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.

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 (Doughty et al., 2014; 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., 2017, 2018; Gaither 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 (Martin et al., 1987; 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 pre-zygotic 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.