Guangpeng Liu
Annalisa Bracco
Andrea M. Quattrini
Santiago Herrera- 1School of Earth and Atmospheric Sciences, Georgia Institute of Technology, Atlanta, GA, United States
- 2Department of Invertebrate Zoology, National Museum of Natural History, Smithsonian Institution, Washington, DC, United States
- 3Department of Biological Sciences, Lehigh University, Bethlehem, PA, United States
Fine-scale larval dispersal and connectivity processes are key to species survival, growth, recovery and adaptation under rapidly changing disturbances. Quantifying both are required to develop any effective management strategy. In the present work, we examine the dispersal pattern and potential connectivity of a common deep-water coral, Paramuricea biscaya, found in the northern Gulf of Mexico by evaluating predictions of physical models with estimates of genetic connectivity. While genetic approaches provide estimates of realized connectivity, they do not provide information on the dispersal process. Physical circulation models can now achieve kilometer-scale resolution sufficient to provide detailed insight into the pathways and scales of larval dispersal. A high-resolution regional ocean circulation model is integrated for 2015 and its advective pathways are compared with the outcome of the genetic connectivity estimates of corals collected at six locations over the continental slope at depths comprised between 1,000 and 3,000 m. Furthermore, the likely interannual variability is extrapolated using ocean hindcasts available for this basin. The general connectivity pattern exhibits a dispersal trend from east to west following 1,000 to 2,000-m isobaths, corresponding to the overall westward near-bottom circulation. The connectivity networks predicted by our model were mostly congruent with the estimated genetic connectivity patterns. Our results show that although dispersal distances of 100 km or less are common, depth differences between tens to a few hundred meters can effectively limit larval dispersal. A probabilistic graphic model suggests that stepping-stone dispersal mediated by intermediate sites provides a likely mechanism for long-distance connectivity between the populations separated by distances of 300 km or greater, such as those found in the DeSoto and Keathley canyons.
Introduction
Deep-water or cold-water corals are long-lived and slow-growing organisms commonly found at depths greater than 50 m (Cairns, 2007; Roark et al., 2009; Sherwood and Edinger, 2009). They play an essential role in providing habitats for a diversity of vertebrate and invertebrate species and are highly susceptible to natural and anthropogenic disturbances (Guinotte et al., 2006; Turley et al., 2007; White et al., 2012; Hoegh-Guldberg et al., 2017). Understanding larval dispersal and connectivity patterns of deep-water corals is a first, necessary step for their conservation and management in response to the multiple threats they face (Palumbi, 2003; Cowen et al., 2007; Botsford et al., 2009).
While shallow corals are generally well-sampled, direct surveys of deep corals are scarce because of the substantial cost and logistical difficulties (Doughty et al., 2014; Quattrini et al., 2015; Girard et al., 2019). The integration of biological data and physical ocean models has helped predict coral habitat suitability (Tong et al., 2013; Hu et al., 2020; Kinlan et al., 2020) and identify larval dispersal and connectivity patterns (Etter and Bower, 2015; Hilario et al., 2015; Breusing et al., 2016; Cardona et al., 2016; Storlazzi et al., 2017; Nolasco et al., 2018; Bracco et al., 2019; Fobert et al., 2019; Gary et al., 2020; Ross et al., 2020). Recent studies suggest that this biophysical framework offers meaningful predictions of connectivity (Gary et al., 2020; Ross et al., 2020), despite the uncertainties related to the sparsity of in situ measurements and model biases in the representation of bottom boundary layer dynamics (see Bracco et al., 2020 for a recent review pertinent to the Gulf of Mexico).
Numerous shallow and deep-water corals populate the northern Gulf of Mexico (GoM) and contribute to the functionality and biodiversity of marine ecosystems (Cordes et al., 2008; Precht et al., 2014; Gil-Agudelo et al., 2020). Deep-water corals in the GoM are subject to various natural and anthropogenic stresses such as increasing water temperatures, acidification, overfishing, and pollution. For example, the 2010 Deep-water Horizon (DWH) oil spill released ∼ 4.1 million barrels of oil into the Gulf (McNutt et al., 2012) and dramatically impacted the vulnerable coral communities in the proximity of the spill site (White et al., 2012; Fisher et al., 2014; Girard et al., 2019). This event and its ecological consequences pointed to the need for restoration actions and improved management strategies, raising interest for a better understanding of the biological and physical processes that affect coral connectivity in the deep-sea.
For deep-water corals, both mesoscale eddies (10–200 Km) and submesoscale circulations (1–10 Km), together with bottom boundary layer turbulence, influence their dispersal (Cardona et al., 2016; Bracco et al., 2019). Near-bottom submesoscale circulations such as fronts, vorticity filaments and small eddies, form due to instabilities induced by shear layers. These circulations can isolate larvae by trapping them inside their cores and transport them to other locations along the continental slope (Bracco et al., 2016). Additionally, the large vertical velocities and diapycnal mixing associated with submesoscale motions can contribute to the vertical transport of the larvae (Bracco et al., 2018; Vic et al., 2018).
Despite the growing number of studies focusing on biophysical dispersal models, connectivity studies that combine genetic data with models resolving the physical circulation and bathymetry at kilometer-scale resolution (submesoscale) are scarce (Cardona et al., 2016; Nolasco et al., 2018; Bracco et al., 2019; Fobert et al., 2019; Gary et al., 2020; Ross et al., 2020). This scarcity is due to the high computational costs of high-spatial resolution models, which limit them short temporal scales (days to months), and the challenges in obtaining sufficiently large sample sizes of deep-sea species for population genetics.
In this work we focus on Paramuricea biscaya, an octocoral in the family Plexauridae. Paramuricea biscaya is one of the most common and abundant corals in the GoM between 1,200 and 2,500 m (Doughty et al., 2014). Populations of this species were directly impacted by the 2010 Deep-water Horizon oil spill (DWH), particularly in the Mississippi Canyon area (White et al., 2012; Fisher et al., 2014), and thus are considered primary targets for restoration (Deepwater Horizon Natural Resource Damage Assessment Trustees, 2016). We investigate the metapopulation connectivity of P. biscaya in the northern GoM using a submesoscale permitting ocean circulation and larval dispersal model. The circulation model resolves the ocean mesoscale circulations (∼10–100 km) and in part the submesoscale (100 m–10 km) range. We also evaluate the performance of the model by comparing potential connectivity probabilities with genetic connectivity estimates.
This paper is a companion to the paper by Galaska et al. (2021) that describes the analyses of genetic connectivity and seascape genomics. Here, we compare the modeled current velocities to those of mesoscale resolving HYCOM-NCODA reanalysis. Furthermore, we explore the factors controlling the larval dispersal pathways and connectivity networks at the sites where P. biscaya occurs, in off-line Lagrangian particle integrations and on-line Eulerian dye simulations through spatial density analysis and a probabilistic graphic model. We also evaluate the potential role of intermediate populations predicted by habitat suitability models (Georgian et al., 2020) as stepping-stones for dispersal. Finally, we discuss the role of annual and inter-annual seasonality in modulating P. biscaya connectivity patterns in the GoM by means of a coarser mesoscale resolving data-assimilative hindcast.
Data and Methods
Large-Scale Circulation of the Study Area
Our study area comprises the northern half of the Gulf of Mexico and includes six sites, named after the lease blocks, where Paramuricea biscaya populations are known. These sites are: DeSoto Canyon 673 (DC673), Mississippi Canyon 344 (MC344), Mississippi Canyon 297 (MC297), Mississippi Canyon 294 (MC294), Green Canyon 852 (GC852), and Keathley Canyon 405 (KC405) (Doughty et al., 2014; Girard et al., 2019; Vohsen et al., 2020; Figure 1). The 2010 Deepwater Horizon oil spill directly impacted P. biscaya populations at the Mississippi Canyon sites MC294, MC297, and MC344 (White et al., 2012; Fisher et al., 2014).
Figure 1. (A) Topography map of the study area showing the sampled and predicted sites hosting P. biscaya populations. (A) Red boxes (0.05 × 0.05°) with colored textbox indicates six main sites (from east to west, DC673, MC344, MC297, MC294, GC852, and KC405). Green boxes with black edge are predicted intermediate suitable habitats selected from models by Georgian et al. (2020). Yellow contours indicate 1,000, 2,000, and 3,000 m isobaths. Light green boxes at the up left corner summary three release strategies in this work. (B) The gray shaded area shows the topographic profile of the northern GoM averaged between 1,000 and 3,000 m. The colored boxes indicate the depth where corals were collected (red) or larval particles released (blue). Gray lines show the average depth of the s-layers in our model.
The large-scale circulation of the study area (shown in Supplementary Figure 1) is dominated in the upper 800–1,000 m by the presence of the anticyclonic Loop Current (LC) that enters the basin through the Yucatan Channel and leaves through the Florida Straits. The LC penetrates northward to about 26.5–27.5°N and usually extends longitudinally east of 86°W (Vukovich, 1988, 2007). Large anticyclonic mesoscale eddies with diameters of about 200 km spin off the main LC at irregular intervals and populate the basin until they dissipate by interacting with the continental slope in the western GoM (Cardona and Bracco, 2016; Donohue et al., 2016). At depth, the large-scale mean circulation is cyclonic (DeHaan and Sturges, 2005; Weatherly et al., 2005). Along the continental slope, bottom currents are highly variable (Supplementary Figure 1), and can intensify due to vortex stretching, submesoscale instabilities and topographic Rossby waves (Hamilton, 2009; Kolodziejczyk et al., 2012; Bracco et al., 2016).
Larval Dispersal Model
Predicting larval dispersal is complicated by the limited knowledge of larvae’s behaviors, especially for deep-sea corals, and by strong, often poorly characterized, variability in deep ocean currents. The application of an integrated biophysical model remains a practicable approach to address this challenge, even though there are notable differences in the estimation of larval travel distance and dispersal pattern among different models (Cowen et al., 2007; Werner et al., 2007; Edmunds et al., 2018; Ross et al., 2020). A typical modeling framework for connectivity studies includes an ocean physical model that provides circulation, or sometimes temperature and salinity, information as background forcing field, and a module for particle (i.e., larvae) tracking those accounts for behavioral characteristics (e.g., age, life span, swimming abilities, larval buoyancy).
In this work, we adopted the three-dimensional Coastal and Regional Ocean Community model (CROCO) that is built upon the Adaptive Grid Refinement in Fortran (AGRIF) version of the Regional Ocean Modeling System (ROMS) (Shchepetkin and McWilliams, 2005; Debreu et al., 2012). It is a split-explicit, hydrostatic, and terrain-following model that is designed for simulating high-resolution nearshore and offshore dynamics and has been used successfully in larval dispersal studies (Kim and Barth, 2011; Cardona et al., 2016; Nolasco et al., 2018; Vic et al., 2018; Bracco et al., 2019; Bani et al., 2020). Here, CROCO covers a large portion of the GoM between 98°–82° W and 24°–31° N (Supplementary Figure 1), and has a grid resolution of about 1 km in the horizontal space and 50 sigma layers in the vertical direction, with the vertical S-coordinate surface stretching parameter set to 5 and the S-coordinate bottom stretching parameter equal to 0.4. Lateral momentum advection is obtained through a 3rd-order upstream biased scheme, and vertical momentum advection through a 4th-order compact advection scheme. The non-linear K-Profile Parameterization (KPP) scheme parameterizes vertical mixing (Large et al., 1994). Three-dimensional tracer advection is achieved through a split and rotated 3rd-order upstream-biased advection scheme in the horizontal, which minimizes spurious diapycnal mixing but does not guarantee positive values of tracer concentration (Marchesiello et al., 2009), and a 4th-order centered advection scheme with harmonic averaging in the vertical. Lateral mixing is parameterized by a Laplacian for momentum and by rotated diffusion with a stabilizing correction for the tracers.
The model bathymetry is derived from the Global Multi-Resolution Topography (GMRT) Synthesis (Ryan et al., 2009) smoothed with a maximum slope factor of 0.25 to reduce horizontal pressure gradient errors (Sikirić et al., 2009). The southern and eastern open boundaries are nudged to the 6-hourly data from the Hybrid Coordinate Ocean Model–Navy Coupled Ocean Data Assimilation (HYCOM-NCODA) Analysis system. Six-hourly atmospheric forcing files (wind stresses, heat fluxes, and daily precipitation) are from the European Centre for Medium-Range Weather Forecast ERA-Interim reanalysis (Poli et al., 2010). Daily freshwater discharges for the five main rivers in the GoM (Mississippi, Atchafalaya, Colorado, Brazos, and Apalachicola) from the United States Geological Survey (USGS) are converted to an equivalent surface freshwater flux that decays away from the river mouths at a constant rate as in Barkan et al. (2017). River momentum flux and tidal forcing are neglected in this work because of their weak influences on the deep-sea area which is the focus of this study (Gouillon et al., 2010; Bracco et al., 2019). Initial conditions are created by interpolating the field of HYCOM on September 31th 2014 to the CROCO grid; the first 4 months of the simulation are discarded as spin-up. CROCO fields are saved every hour for offline particle tracking. At 1 km horizontal resolution, the use of hourly averaged velocity fields introduces only a small error in the tracer advection (Smith et al., 2011; Choi et al., 2017).
Larvae Tracking
Three release experiments with different tracking periods and different types of tracers are conducted in this work to investigate the connectivity pattern of P. biscaya in the northern GoM. Specifically, 4,489 neutrally buoyant Lagrangian tracers (hereafter referred to as larval particles or particles) are deployed uniformly in 0.05 × 0.05° boxes centered at the location of known populations (Supplementary Table 1 and Figure 1A, red box with a colored text box next to it) and at 11 intermediate sites (green boxes) that could host P. biscaya populations according to habitat suitability modeling predictions (Georgian et al., 2020). These particles are tracked off-line (release type 1 and 2) using a Lagrangian tool developed to simulate ichthyoplankton dynamics (Ichthyop) (Lett et al., 2008) and are recorded hourly. Although the actual size of coral larvae is not infinitesimally small and could be slightly negatively buoyant (Miller, 1998; Miller et al., 2012; Brugler et al., 2013), the infinitesimally small approximation holds given the 1 km model resolution. A previous study has shown that in an environment with strong submesoscale features, slightly heavier/lighter (10%) buoyancy does not affect the transport significantly (Zhong et al., 2012). No other biological behaviors such as growth, mortality, settlement, and swimming are considered in this work, given that these are unconstrained for P. biscaya. The CROCO release depths are shown in Figure 1B. There is a 73 m difference on average between observations and model, and the largest discrepancy is found at GC852 (∼200 m), where the observed bathymetry is very steep and varies greatly laterally on scales smaller than the model grid resolution.
A total of 76,313 particles are released in the model layer above that at the seafloor on January 25th, April 25th, July 24th, and November 1st, 2015, and tracked for 56 days (release type 1). The pelagic larval duration (PLD) for Paramuricea biscaya is unknown, but Hilario et al. (2015) found that a PLD between 35 and 69 days may be representative of 50–75% of deep-sea species. In addition, the particles released at the six sampling locations on November 1st are followed for another 92 days to evaluate connectivity over 5 months (∼150 days in total, release type 2). Finally, the evolution of a dye released near the bottom (in the first s-layer) at the six sampling locations is simulated on-line (directly in CROCO) (release type 3) and followed for 120 days to explore the consistency between Lagrangian off-line and Eulerian on-line experiments.
Hybrid Coordinate Ocean Model Hindcast
The inter-annual variability of the near-bottom currents is evaluated using the HYCOM hourly analysis at 1/25° horizontal resolution from 2010 to 2018 (Exp1, HYCOM/GOMl0.04) and the reanalysis data from 2010 to 2012 available at the same horizontal resolution but only at a three-hourly frequency (Exp2, HYCOM/GOMu0.04). The local circulations of these two experiments differ considerably over the common period because of the different choices regarding model configuration, vertical discretization and data assimilation routines (see https://www.hycom.org/hycom/documentation).
We analyzed the velocity field over the whole 2010–2018 period and found that in 2011 the derived near-bottom currents differ significantly between Exp1 and Exp2, and differ the most from those in 2015. Therefore, we chose to simulate the dispersion patterns off-line using the HYCOM data in 2011. Practically, by considering 2015 in CROCO and 2011 in the two HYCOM experiments we are exploring conditions as different as possible within the 2010–2018 period. Particle trajectories in HYCOM are advected using only the 2-dimensional near-bottom velocity field. To make the comparison with our simulations most relevant, we interpolated the original HYCOM velocity field at the horizontal resolution of 5 km and with the same vertical discretization used in the CROCO runs.
Genetic Data
To evaluate the performance of our models in predicting connectivity, we compared our results with the genetic connectivity estimates (migration rates m) by Galaska et al. (2021). Briefly, Galaska et al. (2021) produced single nucleotide polymorphisms (SNPs) data from individuals collected at the six main sites (DC673, MC344, MC297, MC294, GC852, and KC405) using the reduced representation DNA sequencing (RAD-seq) method (Baird et al., 2008; Reitzel et al., 2013). Migration rates (m), defined as the proportion of immigrant individuals in the last two generations, were estimated using BAYESASS v3.0.4.2 (Wilson and Rannala, 2003).
Results and Discussion
Circulation Features in the Gulf of Mexico and Its Annual and Inter-Annual Variability
Figure 2 shows the time-averaged near-bottom zonal velocities over the continental slope between 1,000 and 3,000 m, and their time-series where the coral sites are located, for CROCO in 2015 and for the HYCOM-NCODA analysis from 2010 to 2019. The spatial resolution difference among the two models implies that CROCO partially resolves submesoscale dynamics, while HYCOM-NCODA does not (e.g. Luo et al., 2016). The horizontal patterns of the averaged zonal current (west-east) of HYCOM outputs and our model result in 2015 (Figures 2A,C) illustrate the differences in current speed and, especially, variability (standard deviation) owing to CROCO’s higher resolution (Bracco et al., 2016). At the same time, some similarities are evident. For example, the prevalence of positive/negative velocities around (93.5°W, 27°N)/(93.5°W, 26.5°N), intermittent positive (eastward) and negative (westward) values between 94°W and 90°W, the presence of a westward velocity “belt” following the 3,000 m isobaths. The variability patterns between HYCOM and model outputs are also similar, just of stronger amplitude in CROCO, with the largest variability found between 94.5°W and 93°W, near the Sigsbee Escarpment (89°W–91°W) and the Mississippi Fan (east of 90°W and south of 27°N). Energetic currents and large variability in the vicinity of the Sigsbee Escarpment are supported by field observations (Hamilton and Lugo-Fernandez, 2001). In the De Soto Canyon region (89°W–87°W, 27°N north), on the other hand, currents are weaker and less variable, as documented in previous studies (Bracco et al., 2016; Cardona et al., 2016).
Figure 2. Near-bottom circulation in the northern GoM between 1,000 and 3,000 m. (A) Averaged hourly HYCOM Exp1 zonal current (west–east) from 2010 to 2018, (B) the corresponding standard deviation (std), (C) daily averaged model zonal current, dataset is collected from February, May, August, and November 2015 with a time period of 30 days for each simulation, and (D) the std of the modeled zonal current. The six sampling sites are colored by red dots, while the intermediate sites are shown in green. Panel (E) shows the near-bottom zonal current at each sampling location during the period of 2010–2014 (left), 2015 (middle), and 2016–2018 (right). Gray, black, and color lines indicate result of HYCOM Exp1, HYCOM Exp2, and CROCO, respectively. Positive values indicate an eastward current, and negative values westward.
Figure 2E depicts the time series of zonal near-bottom velocity at the six coral sites. In the period considered, there is some interannual variability but it is not much greater than across different seasons. Currents differ more for amplitude and direction between Exp1 and Exp2, than between Exp1 and CROCO. In 2015, the flow was westward (negative) in the monthly averaged HYCOM output (black thin line) and in the same direction, but generally stronger in CROCO. MC294 is the exception, and the directionality is reversed in HYCOM, but with a small amplitude. Overall in CROCO and Exp1 westward currents prevail around the coral sites.
Mesoscale and submesoscale circulations that may influence the transport of deep-water coral larvae can be seen in the normalized relative vorticity ζ/f plots in Figure 3. ζ = ∂v/∂x –∂u/∂y, where f is the Coriolis parameter. u and v denote the zonal (west-east) and meridional (south-north) velocity components, respectively. x and y are the corresponding distances. The near-bottom vorticity indicates more intense submesoscale circulations over the continental slope in the western side of the domain compared to the eastern one, in agreement with previous work (Cardona et al., 2016). In 2015 slightly stronger submesoscale eddies are found in August compared to the other months, but this is likely the result of interactions between local currents and topography at that time, and is not indicative of a robust seasonal signal. Zoom-in fields of August and November results (Figures 3E,F) provide more details of the interactions between the near-bottom currents and topographic features. Overall, strong cyclonic submesoscale vortices (i.e., positive relative vorticity) are found in “valley” regions outlined by the close seafloor depth contours (see Figure 1A for better visualization) in both seasons. The formation of these small cyclones involves the shear layers and centrifugal instabilities associated with the mean flows and sloping boundary. The generation mechanism is beyond the scope of the present study but is discussed in previous works (Gula et al., 2015; Molemaker et al., 2015; Bracco et al., 2016). Away from these intense cyclonic vortices, GC852 locates in a less stable region with numerous weak, intermittent submesoscale structures.
Figure 3. Relative vorticity near bottom in February (A), May (B), August (C), and November (D) calculated from model simulation. Negative values are anticyclonic or clockwise rotation. Panels (E,F) are the zoom-in regions between 93–90° W, 26.2–27.7° N in August (C) and November (D). Red box with black edge shows the location of GC852.
Genetic Connectivity
The genetic analysis of the Paramuricea samples by Galaska et al. (2021) shows that migration rates are generally low (average m = 0.011), with a few exceptions (Figure 4). Considering the relative locations of sampled sites and the magnitude of the connectivity network, we can infer that the three sites in the Mississippi Canyon (i.e., MC344, 1,853 m; MC297, 1,571 m; and MC294, 1,371 m) are not well-connected to each other (m < 1.5%). This is consistent with the depth-differentiation hypothesis that posits that genetic differentiation is greater across depth than geographic distance, even between sites relatively close by (Quattrini et al., 2015; Bracco et al., 2019; Galaska et al., 2021). Secondly, gene flow is predominantly westward. The DeSoto Canyon site DC673 is a source of genetic material to the Mississippi Canyon sites, particularly the deepest one (MC344). DC673 and KC405 (Keathley Canyon) appear well-connected, despite been separated by 635 km of distance and 522 m of depth. These analyses also indicate that the Green Canyon site GC852 may also be an important source of immigrants to the Mississippi Canyon area because of eastward migration rates.
Figure 4. Connectivity matrix (left) and network (right) calculated from genetic data sampled at the six locations shown in Figure 1A. Figures are modified from Galaska et al. (2021). (A) Matrix values correspond to point estimates of migration rates (m) as percentages (%). (B) Network lines represent connections (dash line for KC405) and dots sites. Dots are color-coded by site. The color of the lines indicates the source site for the connections. Linewidths are proportional to m values. Sites are arranged from east to west (from DC673 to KC405).
Connectivity Pattern in the Physical Circulation Model
The 2015 modeled connectivity patterns among the six sites are illustrated in Figure 5 where the horizontal and vertical distribution of Lagrangian particles is shown after 56 days in each season (Figure 5, release type 1). The horizontal dispersion of particles, even though characterized by detectable differences in the direction of motion and spreading area among sites, is mostly confined between the 1,000 and 2,000 m isobaths in all seasons (Figure 5A). No obvious seasonal dependency is detected except for the relatively wider distribution of particles released from GC852 and KC405 in August, in response to the stronger submesoscale near-bottom flows along the western continental slope compared to the other months (see Figure 3E). In February, May, and November, most particles travel as far as 100–300 km from their release locations, while in August particles released at GC852 and KC405 can be found as far as 600 km.
Figure 5. Distribution of Lagrangian larval particles released at the sampling sites after 56 days. Horizontal (A). Vertical (B). Particles from different sites in panel (A) are colored consistently with previous figures. For each violin in panel (B) the central white square indicates the median, and the bottom and top edges of the black box indicate the interquartile range between 25th and 75th percentiles, respectively. Thin black line shows 95% confidence level, and cyan stars indicate the initial release depth for each site. The width of each violin represents frequency, i.e., density plot.
As mentioned earlier, the near-bottom circulation at the GoM continental slope is predominately along depth contours, therefore larvae migration and connectivity are closely linked to the alongshore (lateral) direction of motion. Virtual larvae released at DC673, MC344, MC297, and MC294 move principally westward in all seasons, as also reported in Cardona et al. (2016). For KC405 and especially GC852 particles, on the other hand, eastward movement can be observed as their lateral velocity is more variable (Figure 2), especially for the February and August releases in the CROCO run. The greater variability of the circulation in the central portion of the GoM continental slope results from the many recirculation zones that occupy this area (see Bracco et al., 2016, their Figures 3, 10).
In the vertical direction (Figure 5B), the particle spreading is highly variable but without a seasonal trend. For the August release, a large portion of particles that originated at GC852 and KC405 are displaced by more than 500 m in 2 months. DC673 is a site characterized by strong eastward and westward currents and high variance in particle displacement (Figures 2A–D) likely caused by the steep slopes of the surrounding topography (Figure 1A). On the contrary, particles released at the Mississippi Canyon sites in all seasons, and GC852 in February, May, and November show smaller variances in their trajectories. For an in-depth discussion of the mechanisms controlling vertical displacement and diapycnal mixing near the ocean bottom in the northern GoM the reader is referred to Bracco et al. (2016, 2020).
The horizontal dispersal patterns for 2015 are quantified by the 2-dimensional kernel density estimation (KDE) using all four releases and the November extended one as well (Figure 6). The KDE is a non-parametric technique to produce a smooth probability density function given a random variable (Rosenblatt, 1956). The figure shows the heat map of particles’ trajectories in the first 10, 30, 56 days and all available releases, i.e., 148 days in November and 56 days in all the other seasons. The last case (Figure 6D) is outside of a biological reasonable scenario, but provides a useful visualization of the connectivity limits in the area of interest. Particles displaced in the vertical direction by more than 800 m depth are discarded (∼20% of all particles and mostly ending in very deep water). For all cases, high KDE values are found in regions within ∼100 km from the release points. A clear pathway from the northeast region near the De Soto Canyon to the southwest area between KC405 and GC852 is outlined following the 2,000 m isobath. In addition to the prevailing westward transport, an eastward branch stemming from KC405 could be responsible for completing the east-to-west connectivity.
Figure 6. Horizontal probability distribution of Lagrangian larval particles after their release at the six sampling sites (colors indicate probability values). Maps of the kernel density estimation (KDE) of particles passages during periods of 10 days (A), 30 days (B), 56 days (C), and all available days (D, 56 days in February, May, and August, and 148 days in November) since release. Subsampling (sample every 20 particles spatially and 2 days temporally) has been applied for plotting. Few larval particles that are displaced more than 800 m in the vertical direction away from the bottom are removed to focus on near-bottom processes.
The potential connectivity network from the model integrations is compared to the observed one from the genetic data in Figure 7. Even though the extended tracking period in November leads to more horizontal (increased from 5 to 6 connections) and vertical (from 2 to 4 connections) connections among the six sampling sites, the modeled network is still less dense than the measured one (genetic), possibly indicating an overall underestimation of coral connectivity. When the vertical dimension is considered, the connectivity from KC405 to GC852 is lost (dash line) and the probability values of several other connections decrease significantly. In line with the genetic evidence, the three MC sites are all connected horizontally, but their connectivity decreases by 64.6% when depth is considered. Connections between the easternmost DC673 and the three MC sites are also observed in the modeled network and are most evident for the DC673-MC344 pair. The model, however, fails to simulate both the long-distance bi-directional communication between DC673 and KC405 and the connectivity out of GC852. In other words, KC405 and GC852 are statistically isolated from the eastern sites in the 2015 CROCO simulation, contrary to the outcome from the genetic inferences.
Figure 7. Horizontal lh (A) and 3D lv (B) connectivity probability (%) matrices and connectivity network (C) calculated from the Lagrangian simulations performed over 56 days in February, May, and August, and 148 days in November. Line colors indicates the source site for each connection. Solid lines indicate three-dimensional connectivity and the dash line represents horizontal connection only. Line width does not correspond to the magnitude of connectivity. Sites are arranged from east to west (from DC673 to KC405).
The Potential Role of Intermediate Populations
The six sampled sites do not represent the only sites hosting P. biscaya populations in the northern GoM. A habitat suitability model recently published by Georgian et al. (2020) predicted a broad distribution of habitat areas where P. biscaya is likely present. Using this model, we selected 11 other potential sites distributed between 1,600 and 2,300 m deep, and roughly equidistantly between the six sampled sites (distances among all sites ranging between approximately 50 and 100 km, Figure 1A). We integrated larvae trajectories from these sites and investigated potential connectivity, under the assumption that these locations are indeed populated by P. biscaya and can thus participate in larval exchange. The connectivity matrices in Figure 8 show the probabilities of the larval exchange among all sites (sampled and predicted). These results clearly show that predicted connectivity is predominantly westward across the region. The probability of connectivity decreases as a function of distance, but depth ultimately dictates whether or not neighboring sites are connected. This is because the connectivity among sites that are relatively close, geographically, is limited by diapycnal mixing across depth.
Figure 8. Heat maps of horizontal lh (A) and 3D lv (B) connectivity probability (%) calculated from all 17 sites (6 sampled and 11 predicted). The connectivity is calculated from 56 days output in February, May, and August and 148 days output in November, 2015. Connectivity values below 0.1% are not shown. Sites are arranged from east to west (from DC673 to KC405).
We built a probabilistic graphic model to quantify the role of intermediate sites in metapopulation connectivity and visualize the modified connectivity matrices. We adopted a directed cyclic graph, instead of the most commonly used Bayesian network, which is a directed acyclic graph (Ben-Gal, 2007). Based on the concepts of conditional probability and chain rule, the joint probability of three events from A (source) to B (intermediate site) to C (sink) can be represented as P(A, B, C) = P(A) × P(B| A) × P(C| A, B). The cyclic network built here allows us to understand the dependency among events (nodes) and assigns probabilities (edges) to them. Our goal was to identify as many connections as possible using the current matrix rather than quantifying the probability/strength of each connection after certain iterations, therefore connections (graph edges) that occur more than once are excluded for the following iterations.
Figure 9 presents the 3-dimensional connectivity pattern and network (based on the vertical matrix in Figure 8) after three iterations/steps (no new connections are found after four or more iterations). Long-distance east to west connections are resolved after three steps (two intermediate sites acting as stepping stones) for particles released from MC344 (Figure 9A). Meanwhile, GC852 particles show both eastward and westward predicted connections despite the direct exchange terminates at an intermediate site around 89°W (Figure 9B). The complete connectivity network (Figure 9C) resolves 19 out of the total 30 possible connections and provides a mechanism for long-distance connectivity between MC344 and KC405 through stepping-stone dispersal mediated by intermediate sites. However, the inferred eastward gene flow from GC852 to the Mississippi Canyon sites remains unexplained by the model.
Figure 9. Predicted 3D patterns of metapopulation connectivity considering the suitable intermediate sites. Connections found after three steps for larval particles released from sites MC344 (A) and GC852 (B). Line colors in panels (A,B) indicate the step at which a connection between two sites was found. Connectivity network among the six sampled sites calculated after including intermediate sites (C). Line colors in panel (C) indicate the source site for each connection. Linewidth in panel (C) indicates direct (thick lines, first step) or mediated connections (regular and thin lines, second and third steps, respectively).
Figure 10. Horizontal kernel density estimation (KDE) pattern of particle locations after 148 days tracking in November without (A) and with (B) extra vertical diffusion (colors indicate probability values). A comparison of the probability density function (PDF) of the depth of particles released at each location is shown in panel (C). Solid/dash line indicates without/with vertical diffusion. All available particles are retained to show a thorough comparison.
Connectivity With Extra Vertical Diffusion
In the work presented so far, vertical and horizontal diffusivities are parameterized at the CROCO grid size. In CROCO, the modeled kz is comparable to observed values (1.3 × 10–4–4 × 10–4 m2 s–1) calculated from a dye injection experiment conducted in 2012 at the Deep-water Horizon site (Ledwell et al., 2016; see e.g., Figure A1 in Bracco et al., 2019). Following previous studies, an additional vertical diffusion coefficient kz of 10–4 m2 s–1 is introduced to the tracking model to recognize the uncertainty associated with larval buoyancy and a likely underestimation of diapycnal mixing very close to the bottom at the continental slope. Since we are exploring the potential impacts of extra perturbation other than simulating the actual marine environment where kz varies spatially (Visser, 1997), a naïve random walk model with constant diffusivity (Hunter et al., 1993) is adopted to track larval particles released in November 2015 for 148 days (Figure 10). Both horizontal KDE and vertical probability density function (PDF) are nearly indistinguishable from the ones obtained without added vertical diffusivity after 5 months. The additional vertical diffusion increases only slightly (by 0.1%) the chances of horizontal connectivity from MC297 to MC344 (not shown).
Lagrangian and Eulerian Tracking Comparison
To further validate the off-line Lagrangian larval particle integrations, we released a dye on-line in the bottom model layer at the six sites in a 0.05 × 0.05° box in November and tracked it for 120 days (release type 3). Figure 11 compares the Lagrangian and Eulerian dispersion based on the distribution of larval particles and the absence/presence of connectivity among sites (Table 1). We stress that the release depths for the Lagrangian larval particles and the dye differ, and that diffusion is included in the dye momentum equations, while this is not the case for the Lagrangian particles (see e.g., Bracco et al., 2009 and more recently Paparella and Vichi, 2020 for pros and cons of using Eulerian vs. Lagrangian approaches). Figure 11 shows the dye concentration field with superimposed particle positions from release type 2 after 120 days of simulation. The results demonstrate that although Lagrangian particles cover a smaller area, they capture the main dispersal features successfully in most cases, e.g., MC344, MC297, GC852, and KC405. For DC673 and MC294, the notable differences are due to the fact that particle positions are plotted at a specific time step (after 120 days), rather than displaying the overall trajectories that would better match the dye pattern. In the end, when comparing connectivity using larval particle trajectories and the 3D dye distribution, it is apparent that the first provides a lower bound to both the horizontal and vertical dispersal due to the lack of subgrid diffusion and not perfectly resolved vertical velocities when using hourly averaged outputs (Wagner et al., 2019).
Figure 11. Horizontal distribution of integrated Eulerian dye concentration normalized by initial tracer concentration (log10 scale, same as Bracco et al., 2018) and Lagrangian larval particles after 120 days from their release on November 1st, 2015 for sampling site at DC673 (A), MC344 (B), MC297 (C), MC294 (D), GC852 (E), and KC405 (F). Colored pentagram with black edge shows release location. Green circles indicate the positions of Lagrangian particles on the same day. Again, particles displaced vertically more than 800 m from the bottom are not included in the calculation to focus on near bottom processes. The number of particles is randomly subsampled for visualization purposes.
Table 1. Difference in Eulerian and Lagrangian connectivity results, H or V indicates that horizontal or vertical connectivity exists in Eulerian but not in Lagrangian result, while HV indicates both horizontal and vertical connectivity are observed in Eulerian but not in Lagrangian result.
Inter-Annual Variability Investigation Using Ocean Hindcast Data
Even in the Eulerian dye approach, the direct connections from GC852 to the Mississippi Canyon site inferred from the genetic data are missing because of the prevailing westward near-bottom current in 2015. While the submesoscale permitting resolution allows for a better representation of the bathymetry and circulation, it increases the computational time and limits the time period of the exploration. To partially address this issue, we examine the Paramuricea connectivity in 2011 when a more persistent eastward circulation was found in the HYCOM hindcast data, especially in Exp2. Using only the near bottom horizontal velocity components for the advection, we released the same number of larval particles as done for CROCO on April 1st 2011, when the mean currents are eastward along the 1,000–2,000 m slopes, and follow them for 148 days. Given the lower frequency at which velocities are saved in Exp2 (3 hourly only), and the underestimation of mean velocities and their standard deviations compared to CROCO results, we also attempted to track the particles by using twice the velocities values reported in the hindcast. Figure 12 provides the horizontal distribution of larval particles after 56-day and 148-day tracking for Exp1 and Exp2 horizontal velocities. The prevalent spreading direction for the particles remained westward, as in CROCO. Greater spreading is found in Exp1 compared to Exp2, and patterns are significantly different. We speculated that the choice of vertical discretization, with better near-bottom vertical resolution in Exp1, explains the differences. Results from Exp1 compare well in terms of overall connectivity patterns with those of CROCO, despite the different times considered (2011 in HYCOM Exp1 and 2015 in CROCO).
Figure 12. Near-bottom distribution of Lagrangian larval particles released in 2-dimensional HYCOM Exp1 (A,B) and HYCOM Exp2 (C,D) fields in April 2011 after 56 days (A,C) and 148 days (B,D) tracking, respectively. Each dot represents the position of a particle. Particles from different sites are colored consistently with previous figures. The setup of vertical layers in each experiment is shown in panel (E).
The traveled distance of the larval particles released at each site in 2015 (CROCO) and 2011 (HYCOM) is shown quantitatively in Figure 13, where the outcome for doubling the horizontal velocities in Exp2 is also displayed. As to be expected, given our selection of a year with conspicuously different currents, inter-annual variability is observed.
Figure 13. Distribution of horizontal dispersal distances for larval particles after 56 and 148 days. November 2015 release using CROCO (A), and April 2011 release using HYCOM for Exp1 (B), Exp2 (C), and Exp2 with doubled zonal velocity (D). Positive values represent eastward movement and negative westward. Particles with a vertical displacement larger than 800 m are removed in CROCO result. The definition of color and violin plot is consistent with Figure 5.
The 75th percentiles of the distance traveled by KC405 particles in the 2011 HYCOM simulation and DC673 particles in the 2015 CROCO run are ∼180 km and ∼300 km (148 days), respectively. These numbers suggest that at least a 10–18 months pelagic larval duration (PLD) is needed to achieve long-distance direct connectivity between KC405 and DC673 if no intermediate sites are considered and without including the vertical aspect. Westward spreading of larval particles in both 2011 and 2015, the generally westward currents in the 9 years considered for HYCOM Exp1, and the weak standard deviation detected in Figure 2 suggest the predominance of an east-to-west along-isobath pathway of dispersal for the three Mississippi Canyon sites. Finally, it is worth noting that by examining the difference between 56- and 148-day results, a fraction of the released particles, e.g., at GC852 (∼17%) in 2015, and at MC344 (∼22–61%) and DC673 (∼70–89%) in 2011, change their traveling directions with respect to their initial locations, indicating a role of seasonal variability.
Summary and Conclusion
In this work, an integrated larval dispersal framework consisting of a high-resolution regional hydrodynamic model (ROMS-CROCO) and a Lagrangian larval particle tracking model (Ichthyop) were performed to predict the dispersal patterns and potential metapopulation connectivity of Paramuricea biscaya in the northern GoM. Lagrangian deployments with 76,313 larval particles were conducted in different seasons for up to ∼150 days and validated by a comparable Eulerian dye experiment in November. The potential contributions of vertical diffusion and intermediate sites to larval connectivity were also studied. The role of inter-annual variability of near-bottom circulations was investigated using HYCOM hindcast data.
The output of our biophysical model showed a mostly congruent agreement with the estimated genetic connectivity patterns (Galaska et al., 2021). In CROCO we found a prevailing westward pathway following the ∼1,000–2,000 m isobath along the continental slope of the northern GoM regardless of seasons in 2015. In general, our estimations of dispersal distances (less than 100 km in 56 days to 300 km in 148 days) agreed well with previous deep-sea studies that considered pelagic larvae duration from 40 days to 1 year (Breusing et al., 2016; Cardona et al., 2016). Strong horizontal but significantly reduced vertical, and therefore 3-dimensional, connectivity among sites near the De Soto Canyon (i.e., DC673, MC344, MC297, and MC294) further confirms the depth differentiation hypothesis in agreement with previous studies (Quattrini et al., 2015; Bracco et al., 2019). The predominantly westward currents and weak variance near the 1,000 m isobath in the De Soto Canyon region found both in HYCOM and in the CROCO model explain the westward confined pathway along the geographic feature shown in Figure 6. In contrast to the relatively stable hydrodynamic environment around the MC sites (also reported in Bracco et al., 2016), highly variable currents over complex topography occupy the central and western portion of the continental slope, around GC852 and KC405 (Figures 1A, 2). Such variability can result at times in eastward transport, as verified in CROCO in February and August 2015 (Figures 5, 6) and in HYCOM in April 2011, with high variance in the modeled displacement of larval particles in both models, and contributes to the diversity of connectivity patterns found in this region.
The inter-annual variability in the near-bottom circulation may be responsible for the few incongruences with the genetic connectivity estimates, for example, the fact that eastward gene flows from GC852 and KC405 is underestimated by the model. Inter-annual variability in the study region has been partially documented by Cardona et al. (2016) in their 3-year simulations from 2010 to 2012. In agreement with our 2011 HYCOM findings, they found several eastward transport events of larval particles originating at ∼92°W to the central Mississippi Canyon (∼89°W) and a predominant westward transport between 600 and 1,000 m depth near the De Soto Canyon.
The long-distance genetic connectivity between DC673 and KC405 may be explained by direct dispersal if we assume a pelagic larval duration of at least 1 year for Paramuricea biscaya. However, this possibility seems unlikely. Pelagic larval durations (PLD) of more than a year have been documented for a few deep-sea invertebrate species (Young et al., 2012), but not for deep-sea corals. Shorter PLDs, between 35 and 69 days, may be representative for most deep-sea species (Hilario et al., 2015). A probabilistic graphic model suggests that stepping-stone dispersal mediated by intermediate sites provides a more likely mechanism for long-distance connectivity between the populations in De Soto and Keathley canyons.
We briefly compared Lagrangian and Eulerian approaches in estimating larval dispersal patterns (Figure 11 and Table 1). Our results implied that the Lagrangian larval particle trajectories computed by interpolating the model’s hourly averaged velocity underestimate vertical velocity and sub-grid diffusion compared to the Eulerian approach (Ali Muttaqi Shah et al., 2017; Wagner et al., 2019). The Lagrangian-derived connectivity represents therefore a lower bound of the Eulerian one. The Lagrangian method presents, on the other hand, several advantages, providing the opportunity to perform multiple sensitivity integrations off-line, and add at low computation cost biological and behavioral constraints, such as mortality, swimming, growing, and settlement. The choice of approach should be done considering the application domain and question(s) in hand.
Our results emphasize the need for multi-year simulations, or at least multi-year analyses of the velocity field, to quantify dispersal patterns in the deep ocean, especially for bi-directional and long-distance connectivity. It is also known that depth and topographic slope are key factors determining the suitability of a habitat for many deep-water corals (Georgian et al., 2020; Hu et al., 2020; Kinlan et al., 2020). Thus, detailed topographic mapping, high horizontal resolution and the fine-scale vertical resolution near the ocean bottom should be adopted to reduce uncertainty in the model representation of bottom currents. Finally, in this work we focused on dispersal processes, but larval traits, e.g., swimming, settlement, and mortality, should be further investigated to improve the realism of modeling studies of coral connectivity.
Data Availability Statement
All numerical tracking outputs, CROCO outputs, and codes for the connectivity analysis used in this research are available at the Dropbox link https://www.dropbox.com/sh/1qrtj8qq9r68ixy/AACacahcanK7mwTX1rTT4laka?dl=0. The Georgia Institute of Technology has a campus-wide license for Dropbox, providing unlimited, ITAR-compliant cloud storage space for large modeling datasets. The connectivity network figures were created using javascript d3.js available at https://d3js.org/, all the other figures were generated with Matlab. CROCO is a community model. The version used in this work is 1.0 and can be found at https://www.croco-ocean.org/download/croco-project/.
Author Contributions
GL and AB conceived the simulations and conducted the analyses. GL performed the model simulations. SH provided the genetic results and analysis. GL, AB, AQ, and SH drafted the manuscript and contributed to the discussions, interpretations, and writing. All authors read and approved the submitted version.
Funding
Funding support for this project was provided by the National Oceanic and Atmospheric Administration’s RESTORE Science Program award NA17NOS4510096 to Lehigh University (SH, AB, and AQ co-PIs). Sampling was supplemented by previous sampling efforts under the Lophelia II Project funded by BOEM and NOAA-OER (BOEM contract #M08PC20038) led by TDI-Brooks International, the NSF RAPID Program (Award #1045079), the NOAA Damage Assessment, Response, and Restoration Program, and ECOGIG (Gulf of Mexico Research Initiative).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We thank Chuck Fisher, Erik Cordes, Illiana Baums for leading those supplemental field efforts and providing access to samples and dive time. We thank Cathy McFadden, Cheryl Morrison, and Frank Parker for project support.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.790927/full#supplementary-material
References
Ali Muttaqi Shah, S. H., Primeau, F. W., Deleersnijder, E., and Heemink, A. W. (2017). Tracing the ventilation pathways of the deep north Pacific Ocean using lagrangian particles and Eulerian tracers. J. Phys. Oceanogr. 47, 1261–1280. doi: 10.1175/jpo-d-16-0098.1
Baird, N. A., Etter, P. D., Atwood, T. S., Currey, M. C., Shver, A. L., Lewis, Z. A., et al. (2008). Rapid SNP discovery and geneticc mapping using sequenced RAD markers. PLoS One 3:e3376. doi: 10.1371/journal.pone.0003376
Bani, R., Marleau, J., Fortin, M.-J., Daigle, R. M., and Guichard, F. (2020). Dynamic larval dispersal can mediate the response of marine metapopulations to multiple climate change impacts. bioRxiv [Preprint] doi: 10.1101/2020.12.05.413062
Barkan, R., McWilliams, J. C., Molemaker, M. J., Choi, J., Srinivasan, K., Shchepetkin, A. F., et al. (2017). Submesoscale dynamics in the northern Gulf of Mexico. Part II: temperature–salinity relations and cross-shelf transport processes. J. Phys. Oceanogr. 47, 2347–2360. doi: 10.1175/jpo-d-17-0040.1
Ben-Gal, I. (2007). Bayesian Networks. Hoboken, NJ: John Wiley & Sons, doi: 10.1002/9780470061572.eqr089
Botsford, L. W., White, J. W., Coffroth, M. A., Paris, C. B., Planes, S., Shearer, T. L., et al. (2009). Connectivity and resilience of coral reef metapopulations in marine protected areas: matching empirical efforts to predictive needs. Coral Reefs 28, 327–337. doi: 10.1007/s00338-009-0466-z
Bracco, A., Choi, J., Joshi, K., Luo, H., and McWilliams, J. C. (2016). Submesoscale currents in the northern Gulf of Mexico: deep phenomena and dispersion over the continental slope. Ocean Model. 101, 43–58. doi: 10.1016/j.ocemod.2016.03.002
Bracco, A., Choi, J., Kurian, J., and Chang, P. (2018). Vertical and horizontal resolution dependency in the model representation of tracer dispersion along the continental slope in the northern Gulf of Mexico. Ocean Model. 122, 13–25. doi: 10.1016/j.ocemod.2017.12.008
Bracco, A., Clayton, S., and Pasquero, C. (2009). Horizontal advection, diffusion, and plankton spectra at the sea surface. J. Geophys. Res 114, 1–4. doi: 10.1029/2007jc004671
Bracco, A., Liu, G., Galaska, M. P., Quattrini, A. M., and Herrera, S. (2019). Integrating physical circulation models and genetic approaches to investigate population connectivity in deep-sea corals. J. Mar. Syst. 198:103189. doi: 10.1016/j.jmarsys.2019.103189
Bracco, A., Paris, C. B., Esbaugh, A. J., Frasier, K., Joye, S. B., Liu, G., et al. (2020). Transport, fate and impacts of the deep plume of petroleum hydrocarbons formed during the macondo blowout. Front. Mar. Sci. 7:5421471. doi: 10.3389/fmars.2020.542147
Breusing, C., Biastoch, A., Drews, A., Metaxas, A., Jollivet, D., Vrijenhoek, R. C., et al. (2016). Biophysical and population genetic models predict the presence of “Phantom” stepping stones connecting Mid-Atlantic Ridge Vent Ecosystems. Curr. Biol. 26, 2257–2267. doi: 10.1016/j.cub.2016.06.062
Brugler, M. R., Opresko, D. M., and France, S. C. (2013). The evolutionary history of the order Antipatharia (Cnidaria: Anthozoa: Hexacorallia) as inferred from mitochondrial and nuclear DNA: implications for black coral taxonomy and systematics. Zool. J. Linnean. Soc. 169, 312–361. doi: 10.1111/zoj.12060
Cairns, S. D. (2007). Deep-water corals: an overview with special reference to diversity and distribution of deep-water scleractinian corals. Bull. Mar. Sci. 81, 311–322.
Cardona, Y., and Bracco, A. (2016). Predictability of mesoscale circulation throughout the water column in the Gulf of Mexico. Deep Sea Res. Part II Top. Stud. Oceanogr. 129, 332–349 doi: 10.1016/j.dsr2.2014.01.008
Cardona, Y., Ruiz-Ramos, D. V., Baums, I. B., and Bracco, A. (2016). Potential connectivity of coldwater black coral communities in the northern Gulf of Mexico. PLoS One 11:e0156257. doi: 10.1371/journal.pone.0156257
Choi, J., Bracco, A., Barkan, R., Shchepetkin, A. F., McWilliams, J. C., and Molemaker, J. M. (2017). Submesoscale dynamics in the northern Gulf of Mexico. Part III: lagrangian implications. J. Phys. Oceanogr. 47, 2361–2376. doi: 10.1175/jpo-d-17-0036.1
Cordes, E. E., McGinley, M. P., Podowski, E. L., Becker, E. L., Lessard-Pilon, S., Viada, S. T., et al. (2008). Coral communities of the deep Gulf of Mexico. Deep Sea Res. Part I Oceanogr. Res. Pap. 55, 777–787. doi: 10.1016/j.dsr.2008.03.005
Cowen, R. K., Gawarkiewicz, G., Pineda, J., Thorrold, S. R., and Werner, F. E. (2007). Population connectivity in marine systems: an overview. Oceanography 20, 14–21. doi: 10.5670/oceanog.2007.2
Debreu, L., Marchesiello, P., Penven, P., and Cambon, G. (2012). Two-way nesting in split-explicit ocean models: algorithms, implem entation and validation. Ocean Model. 49-50, 1–21. doi: 10.1016/j.ocemod.2012.03.003
Deepwater Horizon Natural Resource Damage Assessment Trustees. (2016). Deepwater Horizon Oil Spill: Final Programmatic Damage Assessment and Restoration Plan and Final Programmatic Environmental Impact Statement. Available online at: https://www.gulfspillrestoration.noaa.gov/restoration-planning/gulf-plan (accessed September 5, 2021)
DeHaan, C. J., and Sturges, W. (2005). Deep cyclonic circulation in the Gulf of Mexico. J. Phys. Oceanogr. 35, 1801–1812. doi: 10.1175/JPO2790.1
Donohue, K. A., Watts, D. R., Hamilton, P., Leben, R., and Kennelly, M. (2016). Loop current eddy formation and baroclinic instability. Dyn. Atmos. Oceans 76, 195–216. doi: 10.1016/j.dynatmoce.2016.01.004
Doughty, C. L., Quattrini, A. M., and Cordes, E. E. (2014). Insights into the population dynamics of the deep-sea coral genus Paramuricea in the Gulf of Mexico. Deep Sea Res. Part II Top. Stud. Oceanogr. 99, 71–82. doi: 10.1016/j.dsr2.2013.05.023
Edmunds, P. J., Mcllroy, S. E., Adjeroud, M., Ang, P., Bergman, J., Carpenter, R. C., et al. (2018). Critical information gaps impeding understanding of the role of larval connectivity among coral reef islands in an era of global change. Front. Mar. Sci. 5:290. doi: 10.3389/fmars.2018.00290
Etter, R. J., and Bower, A. S. (2015). Dispersal and population connectivity in the deep North Atlantic estimated from physical transport processes. Deep Sea Res. Part I Oceanogr. Res. Pap. 104, 159–172. doi: 10.1016/j.dsr.2015.06.009
Fisher, C. R., Hsing, P.-Y., Kaiser, C. L., Yoerger, D. R., Roberts, H. H., Shedd, W. W., et al. (2014). Footprint of deep-water horizon blowout impact to deep-water coral communities. Proc. Natl. Acad. Sci. U.S.A. 111, 11744–11749. doi: 10.1073/pnas.1403492111
Fobert, E. K., Treml, E. A., and Swearer, S. E. (2019). Dispersal and population connectivity are phenotype dependent in a marine metapopulation. Proc. Biol. Sci. 286:20191104. doi: 10.1098/rspb.2019.1104
Galaska, M., Liu, G., West, D., Erickson, K., Quattrini, A. M., Bracco, A., et al. (2021). Seascape genomics reveals metapopulation connectivity network of Paramuricea biscaya in the northern Gulf of Mexico. Front Marine Sci. 8:790929. doi: 10.1101/2021.10.06.463359
Gary, S. F., Fox, A. D., Biastoch, A., Roberts, J. M., and Cunningham, S. A. (2020). Larval behaviour, dispersal and population connectivity in the deep sea. Sci. Rep. 10:10675. doi: 10.1038/s41598-020-67503-7
Georgian, S. E., Kramer, K., Saunders, M., Shedd, W., Roberts, H., Lewis, C., et al. (2020). Habitat suitability modelling to predict the spatial distribution of cold-water coral communities affected by the deep-water horizon oil spill. J. Biogeogr. 47, 1455–1466. doi: 10.1111/jbi.13844
Gil-Agudelo, D. L., Cintra-Buenrostro, C. E., Brenner, J., González-Díaz, P., Kiene, W., Lustic, C., et al. (2020). Coral reefs in the Gulf of Mexico large marine ecosystem: conservation status, challenges, and opportunities. Front. Mar. Sci. 6:807. doi: 10.3389/fmars.2019.00807
Girard, F., Cruz, R., Glickman, O., Harpster, T., and Fisher, C. R. (2019). In situ growth of deep-sea octocorals after the deep-water horizon oil spill. Elem. Sci. Anth. 7:19. doi: 10.1525/journal.elementa.349
Gouillon, F., Morey, S. L., Dukhovskoy, D. S., and O’Brien, J. J. (2010). Forced tidal response in the Gulf of Mexico. J. Geophys. Res. Oceans 115:C01002. doi: 10.1029/2010JC006122
Guinotte, J. M., Orr, J., Cairns, S., Freiwald, A., Morgan, L., and George, R. (2006). Will human-induced changes in seawater chemistry alter the distribution of deep-sea scleractinian corals? Front. Ecol. Environ. 4, 141–146. doi: 10.1890/1540-92952006004[0141:WHCISC]2.0.CO;2
Gula, J., Molemaker, M. J., and McWilliams, J. C. (2015). Topographic vorticity generation, submesoscale instability and vortex street formation in the Gulf stream. Geophys. Res. Lett. 42, 4054–4062. doi: 10.1002/2015gl063731
Hamilton, P. (2009). Topographic rossby waves in the Gulf of Mexico. Prog. Oceanogr. 82, 1–31. doi: 10.1016/j.pocean.2009.04.019
Hamilton, P., and Lugo-Fernandez, A. (2001). Observations of high speed deep currents in the northern Gulf of Mexico. Geophys. Res. Lett. 28, 2867–2870. doi: 10.1029/2001gl013039
Hilario, A., Metaxas, A., Gaudron, S. M., Howell, K. L., Mercier, A., Mestre, N. C., et al. (2015). Estimating dispersal distance in the deep sea: challenges and applications to marine reserves. Front. Mar. Sci. 2:6. doi: 10.3389/fmars.2015.00006
Hoegh-Guldberg, O., Poloczanska, E. S., Skirving, W., and Dove, S. (2017). Coral reef ecosystems under climate change and ocean acidification. Front. Mar. Sci 4:158. doi: 10.3389/fmars.2017.00158
Hu, Z., Hu, J., Hu, H., and Zhou, Y. (2020). Predictive habitat suitability modeling of deepsea framework-forming scleractinian corals in the Gulf of Mexico. Sci. Total Environ. 742:140562. doi: 10.1016/j.scitotenv.2020.140562
Hunter, J. R., Craig, P. D., and Phillips, H. E. (1993). On the use of random walk models with spatially variable diffusivity. J. Comput. Phys. 106, 366–376. doi: 10.1016/S0021-9991(83)71114-9
Kim, S., and Barth, J. A. (2011). Connectivity and larval dispersal along the oregon coast estimated by numerical simulations. J. Geophys. Res. 116:C06002. doi: 10.1029/2010jc006741
Kinlan, B. P., Poti, M., Drohan, A. F., Packer, D. B., Dorfman, D. S., and Nizinski, M. S. (2020). Predictive modeling of suitable habitat for deep-sea corals offshore the Northeast United States. Deep Sea Res. Part I Res. Pap. 158:103229. doi: 10.1016/j.dsr.2020.103229
Kolodziejczyk, N., Ochoa, J., Candela, J., and Sheinbaum, J. (2012). Observations of intermittent deep currents and eddies in the Gulf of Mexico. J. Geophys. Res. Oceans 117:9014 doi: 10.1029/2012jc007890
Large, W. G., McWilliams, J. C., and Doney, S. C. (1994). Oceanic vertical mixing: a review and a model with a nonlocal boundary layer parameterization. Rev. Geophys 32, 363–403. doi: 10.1029/94rg01872
Ledwell, J. R., He, R., Xue, Z., DiMarco, S. F., Spencer, L. J., and Chapman, P. (2016). Dispersion of a tracer in the deep Gulf of Mexico. J. Geophys. Res. Oceans 121, 1110–1132. doi: 10.1002/2015JC011405
Lett, C., Verley, P., Mullon, C., Parada, C., Brochier, T., Penven, P., et al. (2008). A Lagrangian tool for modelling ichthyoplankton dynamics. Environ. Model. Softw. 23, 1210–1214. doi: 10.1016/j.envsoft.2008.02.005
Luo, H., Bracco, A., Cardona, Y., and McWilliams, J. C. (2016). Submesoscale circulation in the northern Gulf of Mexico: surface processes and the impact of the freshwater river input. Ocean Model. 101, 68–82. doi: 10.1016/j.ocemod.2016.03.003
Marchesiello, P., Debreu, L., and Couvelard, X. (2009). Spurious diapycnal mixing in terrain-following coordinate models: the problem and a solution. Ocean Model. 26, 156–169. doi: 10.1016/j.ocemod.2008.09.004
McNutt, M. K., Camilli, R., Crone, T. J., Guthrie, G. D., Hsieh, P. A., Ryerson, T. B., et al. (2012). Review of flow rate estimates of the deep-water horizon oil spill. Proc. Natl. Acad. Sci. U.S.A. 109, 20260–20267. doi: 10.1073/pnas.1112139108
Miller, K. J. (1998). Short-distance dispersal of black coral larvae:inference from spatial analysis of colony genotypes. Mar. Ecol. Prog. Ser. 163, 225–233. doi: 10.3354/meps163225
Miller, R. J., Hocevar, J., Stone, R. P., and Fedorov, D. V. (2012). Structure-forming corals and sponges and their use as fish habitat in Bering Sea submarine canyons. PLoS One 7:e33885. doi: 10.1371/journal.pone.0033885
Molemaker, M. J., McWilliams, J. C., and Dewar, W. K. (2015). Submesoscale instability and generation of mesoscale anticyclones near a separation of the California undercurrent. J. Phys. Oceanogr. 45, 613–629. doi: 10.1175/jpo-d-13-0225.1
Nolasco, R., Gomes, I., Peteiro, L., Albuquerque, R., Luna, T., Dubert, J., et al. (2018). Independent estimates of marine population connectivity are more concordant when accounting for uncertainties in larval origins. Sci. Rep. 8:2641. doi: 10.1038/s41598-018-19833-w
Palumbi, S. R. (2003). Population genetics, demographic connectivity, and the design of marine reserves. Ecol. Appl 13, 146–158. doi: 10.1890/1051-07612003013[0146:PGDCAT]2.0.CO;2
Paparella, F., and Vichi, M. (2020). Stirring, mixing, growing: microscale processes change larger scale phytoplankton dynamics. Front. Mar. Sci. 7:654. doi: 10.3389/fmars.2020.00654
Poli, P., Healy, S. B., and Dee, D. P. (2010). Assimilation of global positioning system radio occultation data in the ECMWF ERA-interim reanalysis. Q. J. R. Meteorolog. Soc. 136, 1972–1990. doi: 10.1002/qj.722
Precht, W. F., Deslarzes, K. J. P., Hickerson, E. L., Schmahl, G. P., Nuttall, M. F., and Aronson, R. B. (2014). Back to the future: the history of acroporid corals at the flower garden banks, Gulf of Mexico, USA. Mar. Geol. 349, 152–161. doi: 10.1016/j.margeo.2013.12.012
Quattrini, A. M., Baums, I. B., Shank, T. M., Morrison, C. L., and Cordes, E. E. (2015). Testing the depth-differentiation hypothesis in a deep-water octocoral. Proc. Biol. Sci. 282:20150008. doi: 10.1098/rspb.2015.0008
Reitzel, A. M., Herrera, S., Layden, M. J., Martindale, M. O., and Shank, T. M. (2013). Going where traditional markers have not gone before: utility of and promise for RAD sequencing in marine invertebrate phylogeography and population genomics. Mol. Ecol. 22, 2953–2970. doi: 10.1111/mec.12228
Roark, E. B., Guilderson, T. P., Dunbar, R. B., Fallon, S. J., and Mucciarone, D. A. (2009). Extreme longevity in proteinaceous deep-sea corals. Proc. Natl. Acad. Sci. U.S.A. 106, 5204–5208.
Rosenblatt, M. (1956). Remarks on some nonparametric estimates of a density function. Ann. Math. Stat. 27, 832–837. doi: 10.1214/aoms/1177728190
Ross, R. E., Nimmo-Smith, W. A. M., Torres, R., and Howell, K. L. (2020). Comparing deep-sea larval dispersal models: a cautionary tale for ecology and conservation. Front. Mar. Sci. 7:431. doi: 10.3389/fmars.2020.00431
Ryan, W. B. F., Carbotte, S. M., Coplan, J. O., O’Hara, S., Melkonian, A., Arko, R., et al. (2009). Global multi-resolution topography synthesis. Geochem. Geophys. Geosyst. 10:Q03014. doi: 10.1029/2008gc002332
Shchepetkin, A. F., and McWilliams, J. C. (2005). The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean Model. 9, 347–404. doi: 10.1016/j.ocemod.2004.08.002
Sherwood, O. A., and Edinger, E. N. (2009). Ages and growth rates of some deep-sea gorgonian and antipatharian corals of Newfoundland and Labrador. Can. J. Fish. Aquat.Sci. 66, 142–152. doi: 10.1139/f08-195
Sikirić, M. D., Janekoviæ, I., and Kuzmiæ, M. (2009). A new approach to bathymetry smoothing in sigma-coordinate ocean models. Ocean Model. 29, 128–136. doi: 10.1016/j.ocemod.2009.03.009
Smith, K. S., Keating, S. R., and Kramer, P. R. (2011). Diagnosing lateral mixing in the upper ocean with virtual tracers: spatial and temporal resolution dependence. J. Phys. Oceanogr. 41, 1512–1534. doi: 10.1175/2011jpo4580.1
Storlazzi, C. D., van Ormondt, M., Chen, Y.-L., and Elias, E. P. L. (2017). Modeling fine-scale coral larval dispersal and interisland connectivity to help designate mutually-supporting coral reef marine protected areas: insights from Maui Nui, Hawaii. Front. Mar. Sci. 4:381. doi: 10.3389/fmars.2017.00381
Tong, R., Purser, A., Guinan, J., and Unnithan, V. (2013). Modeling the habitat suitability for deep-water gorgonian corals based on terrain variables. Ecol. Inf. 13, 123–132. doi: 10.1016/j.ecoinf.2012.07.002
Turley, C. M., Roberts, J. M., and Guinotte, J. M. (2007). Corals in deep-water: will the unseen hand of ocean acidification destroy cold-water ecosystems? Coral Reefs 26, 445–448. doi: 10.1007/s00338-007-0247-5
Vic, C., Gula, J., Roullet, G., and Pradillon, F. (2018). Dispersion of deep-sea hydrothermal vent effluents and larvae by submesoscale and tidal currents. Deep Sea Res. Part I Oceanogr. Res. Pap. 133, 1–18. doi: 10.1016/j.dsr.2018.01.001
Visser, A. W. (1997). Using random walk models to simulate the vertical distribution of particles in a turbulent water column. Mar. Ecol. Prog. Ser. 158, 275–281.
Vohsen, S. A., Gruber-Vodicka, H. R., Osman, E. O., Saxton, M. A., Joye, S. B., Dubilier, N., et al. (2020). Deep-sea corals near cold seeps associate with chemoautotrophic bacteria that are related to the symbionts of cold seep and hydrothermal vent mussels. bioRxiv [Preprint]. doi: 10.1101/2020.02.27.968453
Vukovich, F. M. (1988). Loop current boundary variations. J. Geophys. Res. 93:15585. doi: 10.1029/JC093iC12p15585
Vukovich, F. M. (2007). Climatology of ocean features in the Gulf of Mexico using satellite remote sensing data. J. Phys. Oceanogr. 37, 689–707. doi: 10.1175/jpo2989.1
Wagner, P., Rühs, S., Schwarzkopf, F. U., Koszalka, I. M., and Biastoch, A. (2019). Can lagrangian tracking simulate tracer spreading in a high-resolution ocean general circulation model? J. Phys. Oceanogr. 49, 1141–1157. doi: 10.1175/jpo-d-18-0152.1
Weatherly, G., Wienders, N., and Romanou, A. (2005). Intermediate-depth circulation in the Gulf of Mexico estimated from direct measurements, Washington DC. Am. Geophys. Union Geophys. Monogr. Ser. 161, 315–324. doi: 10.1029/161GM22
Werner, F. E., Cowen, R. K., and Paris, C. B. (2007). Coupled biological and physical models: present capabilities and necessary developments for future studies of population connectivity. Oceanography 20, 54–69. doi: 10.5670/oceanog.2007.29
White, H. K., Hsing, P.-Y., Cho, W., Shank, T. M., Cordes, E. E., Quattrini, A. M., et al. (2012). Impact of the deepwater horizon oil spill on a deep-water coral community in the Gulf of Mexico. Proc. Natl. Acad. Sci. U.S.A. 109, 20303–20308. doi: 10.1073/pnas.1118029109
Wilson, G. A., and Rannala, B. (2003). Bayesian inference of recent migration rates using multilocus genotypes. Genetics 163, 1177–1191.
Young, C. M., He, R., Emlet, R. B., Li, Y., Qian, H., Arellano, S. M., et al. (2012). Dispersal of deep-sea larvae from the intra-American seas: simulations of trajectories using ocean models. Integr. Comp. Biol. 52, 483–496.
Keywords: genetic analysis, Gulf of Mexico, connectivity, dispersal models, deep-sea
Citation: Liu G, Bracco A, Quattrini AM and Herrera S (2021) Kilometer-Scale Larval Dispersal Processes Predict Metapopulation Connectivity Pathways for Paramuricea biscaya in the Northern Gulf of Mexico. Front. Mar. Sci. 8:790927. doi: 10.3389/fmars.2021.790927
Received: 07 October 2021; Accepted: 17 November 2021;
Published: 10 December 2021.
Edited by:
Yoshihiro Fujiwara, Japan Agency for Marine-Earth Science and Technology (JAMSTEC), JapanReviewed by:
Jeff Eble, Florida Institute of Technology, United StatesJesus Dubert, University of Aveiro, Portugal
Copyright © 2021 Liu, Bracco, Quattrini and Herrera. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Guangpeng Liu, gliu87@gatech.edu