Exploitation Drives Changes in the Population Connectivity of Queen Conch (Aliger gigas)

The queen conch, Aliger gigas, is an endemic and iconic marine gastropod of the Wider Caribbean region that has been harvested for thousands of years. Conch are slow-moving and require contact to mate; overfishing has reduced populations in many areas compromising its rates of reproduction. Long-range dispersal and mixing between distinct populations occur in the queen conch’s early life history stages, when pelagic larvae are transported by oceanic currents. Genetic studies suggest that gene flow between populations decreases as the distance between populations increases. Here, we assessed how the population connectivity of conch changes with spatially variable patterns of fishing exploitation by simulating larval dispersal and comparing the potential connectivity under an unexploited and a contemporary exploited reproductive scenario. Results demonstrate that reduced egg production, due to heterogeneous fishing pressure and localized depletion, significantly alters population connectivity patterns as well as the structuring of populations and metapopulations across the species’ range. This strongly suggests that estimates of contemporary demographic rates, together with estimates of reproductive output need to be included in population connectivity studies. The existence of self-sustained metapopulations of queen conch throughout the Wider Caribbean suggests that replenishment through larval dispersal occurs primarily within sub-regional spatial scales, emphasizing the need for regional and local conservation and management measures to build and protect reproductively active populations and nursery habitat across multiple jurisdictions.


INTRODUCTION
The queen conch (Aliger gigas, Maxwell et al., 2020) is an iconic marine gastropod endemic to the Wider Caribbean Region (WCR).The species has been harvested for subsistence, economic, and cultural purposes for thousands of years.The queen conch has been listed in Appendix II of the Convention on International Trade in Endangered Species (CITES) since 1992 due to the significant decline of its populations throughout the Caribbean.Intense fisheries and illegal harvest of queen conch have had a significant impact on population recovery (Stoner and Appeldoorn, 2021).Management of conch fisheries is highly variable across the Caribbean region, including strategies such as marine protected areas, minimum shell lengths and/or lip thicknesses for harvest, fisheries quotas, gear restrictions, and closed seasons (Prada et al., 2017).However, despite these regulations, populations appear to be overfished and in a continuing decline in many Caribbean locations (Stoner and Appeldoorn, 2021).
Adult conch are found over a wide range of habitats, particularly on sand flats and grass beds, but also on rubble and gravel, mostly at depths of 20 meters or shallower (Stoner et al., 1994;McCarthy, 2007;Stoner and Davis, 2010).Queen conch, as all gastropods, are slow-moving organisms (Doerr and Hill, 2013) and require direct physical contact for reproduction (Randall, 1964).Consequently, adult (i.e., mature) density is an important factor governing mating success and reproductive output.The principle of inverse density dependence at low densities is referred to as an Allee effect (Allee, 1931).Recent simulation models (Farmer and Doerr, 2022) indicate that the lack of mating observed at low densities may be attributed to low encounter rates of adults (Stoner and Ray-Culp, 2000) and limited sexual facilitation leading to delayed functional maturity (Gascoigne and Lipcius, 2004).When reproductive fitness declines such that the per capita population growth rate becomes negative, localized extinction may result (Allee, 1931;Courchamp et al., 1999).Longterm observations of depleted populations indicate that even with total fishing moratoria, recovery may take decades or longer (Stoner and Appeldoorn, 2021).Consequently, it is critical for fisheries management to identify how localized depletion patterns, which decrease conch densities and reproductive output, may affect the population throughout its entire distribution.
Unfortunately, inconsistent fishery-independent surveys using different methodologies only provide limited information on current conch densities throughout the WCR (Prada et al., 2017).While there is genetic evidence of stock mixing throughout the Caribbean with some localized population isolation patterns (Peŕez-Enriquez et al., 2011;Zamora-Bustillos et al., 2011;Maŕquez-Pretel et al., 2013;Machkour-M'Rabet et al., 2017;Truelove et al., 2017;Blythe-Mallett et al., 2021), little is known about the ecological connectivity of queen conch throughout its range.The species is largely sessile; thus, most long-distance connections occur during the planktonic early life stages (Davis, 2005;Kough et al., 2019).Indeed, spatial genetic structure has been detected at different spatial scales, both between and within territorial waters of Caribbean jurisdictions, suggesting that queen conch do not form a single panmictic population in the greater Caribbean and that variation in marine currents, surface winds, and meteorological events can either promote larval dispersal or act as barriers enhancing local larval retention (Peŕez-Enriquez et al., 2011;Zamora-Bustillos et al., 2011;Maŕquez-Pretel et al., 2013;Machkour-M'Rabet et al., 2017;Truelove et al., 2017;Kitson-Walters et al., 2018;Blythe-Mallett et al., 2021).While modeling studies of larval dispersal of queen conch corroborate these findings at regional scales (Kough et al., 2019;Souza and Kough, 2020), these studies focused on elucidating connectivity patterns mostly around The Bahamas.Thus, it is still unknown how the interplay of physical processes, biological traits, and localized depletion patterns affect the demographic patterns of queen conch throughout the WCR.Moreover, it is unknown if contemporary adult densities, emerging particularly from localized exploitation patterns, can affect population-wide demographics, as suggested for some populations (Zamora-Bustillos et al., 2011;Maŕquez-Pretel et al., 2013).
Here we aim to understand how fisheries can alter the population structure and demographics of the queen conch by estimating ecological connectivity under an "unexploited" historical density scenario and an "exploited" contemporary scenario that incorporates realistic localized density patterns.Estimated connectivity patterns are compared among different scenarios of reproductive output to evaluate how localized harvest patterns can affect the population-wide demographics of the species.This is the first effort to identify connectivity patterns of the species throughout the Wider Caribbean Region and incorporate viable reproduction using a comprehensive set of adult densities.The results should inform collaborative management among the various Caribbean states and territories for one of the region's most economically important species.

Connectivity Modeling System
To estimate the population connectivity of queen conch, we used the Connectivity Modeling System (CMS), a community multiscale modeling system, with an Individual-Based Model and stochastic Lagrangian framework (Paris et al., 2013).The CMS simulates complex biophysical dispersive processes during the pelagic early life history of marine organisms, producing estimates of probabilistic connectivity among populations.Our simulations thus captured the queen conch's early pelagic life-history stages, from hatching to the settlement of competent veligers.
In CMS, embryos and larvae are advected by deterministic velocity fields from hydrodynamic models and by a stochastic component representing unresolved sub-grid scale processes.The velocity field is interpolated for each individual location and time step using a tricubic scheme and integrated to the next time step with a unique spatiotemporal forth-order Runge-Kutta scheme; the stochastic movement is given by a random walk model, scaled by the diffusivity coefficient of the system (Okubo, 1971).CMS also simulates biological traits affecting an organism's dispersal, such as ontogenetic vertical migration (OVM).Spawning and settlement habitats are represented in the model by 2D discrete polygons.Larvae are advected for the time defined as their maximum pelagic larval duration (PLD), and those older than a defined competency period that are found within polygon boundaries are considered settled and removed from the simulation.The CMS saves individual larval positions and attributes through time, together with files with the spawning and settlement sites and times of successful settlers, allowing for the estimation of dispersal and connectivity patterns.

Ocean Velocity Fields
We used velocity fields from two well-established hydrodynamic products, the hindcast and reanalysis of the Global HYCOM and the Mercator GLORYS12v1 reanalysis (Figure 1, hereafter HYCOM and GLORYS, respectively).By considering different hydrodynamic models, we could capture stochastic changes in flow patterns and explicitly address the uncertainty of connectivity estimates due to the underlying hydrodynamic flow fields.Our model domain extended from 8°N to 37°N and from 98°W to 59°W , and our temporal range was from March to December of 2013 to 2017, capturing both the spatial distribution of the conch and the seasonality of its spawning season.The five-year range is adequate to capture interannual variability in the physical oceanography and stochastic patterns of dispersal.
The Mercator GLORYS version 12v1 is the reanalysis from the Copernicus Marine Environment Monitoring Service (CMEMS), largely based on the CMEMS real-time forecast system (Lellouche et al., 2018;Lellouche et al., 2021).It presents global coverage at 1/12°(ca.8km) horizontal resolution and 50 vertical levels, and it is available at daily intervals for our study period.The model component is version 3.1 of the NEMO (Nucleus for European Modelling of the Ocean) model, driven at surface by ECMWF ERA-Interim reanalysis.It uses a reduced-order Kalman filter to assimilate sea level anomaly, sea surface temperature, and sea ice concentration from satellite measurements, and in situ temperature and salinity vertical profiles from the CORA (Coriolis Ocean Database ReAnalysis, Lellouche et al., 2018).It also corrects large-scale bias in temperature and salinity.

Habitat
The majority of queen conch spawning habitat occurs in waters less than 20 m depth, in areas of coarse sand and rubble, sometimes with mixed seagrass (Salley, 1986;Berg Jr. et al., 1992;Boidron-Metairon, 1992;Stoner and Sandt, 1992;Stoner and Schwarte, 1994;Delgado and Glazer, 2020).Nursery habitat for queen conch is found in shallower sheltered back reef areas, primarily associated with native seagrass meadows (Stoner, 1989;Stoner et al., 1994;Stoner et al., 1996;Jones and Stoner, 1997;Doerr and Hill, 2018;Boman et al., 2019).Juveniles are mostly found in seagrass but can also recruit to shallow habitats such as algal flats (Glazer and Berg Jr., 1994).
We compared our habitat area estimates to published seagrass habitat cover and conch fishing areas (SI Figure 3) including compilations of global geomorphic zones (UNEP-WCMC and Short 2021; Allen Coral Atlas, 2020;McKenzie et al., 2020;Schill et al., 2021); studies focused on jurisdictions or regional levels (Wabnitz et al., 2008;Tewfik et al., 2017;Leoń-Peŕez et al., 2019) and fishing sites (compiled in Prada et al., 2017).Estimates of seagrass area by jurisdiction were highly variable, and estimates of conch fishing areas were generally much lower than the highest estimates of seagrass cover.Overall, our habitat area estimates were much lower than total seagrass area estimates per jurisdiction and ranged from ~30% to 100% of the area of the fishing bank estimates (SI Figure 3).Our habitat areas were thus a conservative measurement of conch habitat throughout the Caribbean.

Adult Densities
Data for adult densities were obtained from a recent comprehensive literature review which compiled fishery-independent densities of adult conch populations for each of the 42 jurisdictions in our study area (Horn et al., 2022, Figure 3).Jurisdiction here was defined as any state, territory, collectivity, municipality, or country within the species range having management authority over queen conch stocks.For each of the 42 jurisdictions, all available comprehensive survey-level estimates of adult queen conch per hectare were compiled to generate cross-shelf adult conch density estimates.
Data availability was variable among jurisdictions, and some estimates reflect recent studies which conducted randomized FIGURE 3 | Estimated adult conch densities for the jurisdictions considered in this study (reproduced from Horn et al. 2022, sources listed in Appendix 1).Data points are sized relative to densities; green symbols indicate conch populations with >100 adult conch/ha, gold symbols indicate 50-99.9adult conch/ha, and red symbols indicate <50 adult conch/ha.Countries without density data are indicated with an X -data for these countries were interpolated from the nearest neighbor.Gray circles represent the 95% confidence interval (2.5 -97.5% bounds) of density estimates for each jurisdiction.Where densities are summarized on a jurisdiction level, the points appear at the approximate center point of the jurisdiction; where densities are summarized on a subregional level the points appear at the location of fishing banks or subregions (i.e., The Bahamas and Puerto Rico).Please refer to Horn et al. 2022 for more details on the methodology.
surveys over multiple habitat and shelf areas, while others reflect punctual surveys or data from workshops (SI Appendix 1, Horn et al., 2022).Cross-shelf surveys likely generate more suitable estimates of overall queen conch populations, and cross-shelf surveys seem to be a widely used monitoring method for queen conch stocks.Thus, intra-aggregation density data (e.g., from the few surveys conducted within spawning aggregations) were removed from the dataset (SI Appendix 1).
Priority was given to using surveys carried out after 2000; when this was not possible (Bermuda, Dominican Republic), surveys from the 1990s were used.For jurisdictions where no data were available, densities from the nearest neighbor jurisdiction were used (Costa Rica, Curacao, Dominica, Grenada, Montserrat, Saint Kitts and Nevis, Saint-Barthelemy, Trinidad and Tobago).When available, we also considered available information regarding the status of stocks and local exploitation/consumption of conch.However, for most jurisdictions, there were multiple independent estimates available (SI Appendix 1).In such cases, we fitted a distribution to the independent estimates of mean conch density (see Horn et al., 2022 for details) and then used the median of the distribution as a best estimate of mean adult queen conch density for each jurisdiction.For jurisdictions with large shelf areas and with data available across multiple regions, significant differences across regions were assessed.Densities for The Bahamas, Colombia, and Puerto Rico were significantly different among regions, and separate density estimates were created for these areas.For The Bahamas, the difference is considered between regional banks (Little Bahama Bank, Western Great Bahama Bank, Central Great Bahama Bank, Jumentos and Ragged, and Cay Sal); for Colombia, the division is between nearshore sites and offshore banks located hundreds of miles away (Serranilla, Quitasueño and Serrana Banks); and for Puerto Rico is between shallow sites and a putative mesophotic spawning site at the Abrir La Sierra mesophotic reef (located in the southeast of the Mona Passage, Garcia-Sais et al., 2012).
The compiled density (Figure 3, Horn et al., 2022) showed that most jurisdictions, or regional divisions of jurisdictions (31 out of 46 jurisdictions), present populations with adult densities varying from 0-50 conch/ha, four jurisdictions present densities from 50-100 conch/ha, and 11 present adult densities over 100 conch/ha.

Spawning Behavior
Queen conch exhibit a seasonal spawning pattern, which varies with latitude across our study region.In our model approach we used a standardized annual spawning curve for all jurisdictions, varying the month of onset and peak of spawning for each jurisdiction.Annual seasonal cycles of spawning activity were derived from five studies which evaluated annual egg production (Weil and Laughlin, 1984;Stoner et al., 1992;Appeldoorn, 1993;Delgado and Glazer, 2020), and each was fit to a series of curves (Normal, Weibull, Lognormal, and Gamma).The Weibull curve was the best fit overall, and the best-fit Weibull curve parameters were then averaged to create a single curve that could be applied to all jurisdictions, representing that the ramping up of spawning occurs at a slightly higher rate than the ramping down of spawning after the peak (SI Figure 1).We compiled the month of spawning peak and onset for each jurisdiction from published visual surveys of reproductive output and gonadal histological studies (Randall, 1964;D'Asaro, 1965;Hesse, 1976;Brownell, 1977;Davis et al., 1984;Appeldoorn, 1985;Cruz, 1986;Salley, 1986;Corral and Ogawa, 1987;Coulston et al., 1987;Berg et al., 1992;Stoner et al., 1992, Appeldoorn, 1993. Márquez-Pretel et al., 1994;Aldana-Aranda et al., 2003;Avila-Poveda and Baqueiro-Caŕdenas, 2009;Bissada, 2012;Appeldoorn and Baker, 2013;Cala et al., 2013;Aldana-Aranda et al., 2014;Meijer, 2014;de Graaf et al., 2015;Wynee, 2016;Boman et al. 2018;Tiley et al. 2018;Delgado and Glazer, 2020).Expected spawning activity was then interpolated over all months of the year by using the Weibull curve representing the annual cycle, shifting it temporally such that the peak of the curve matched the reported peak.The expected spawning activity was thus expressed as a series of percentages by month, with the sum of all months equal to 100%.This extrapolation method produced patterns of spawning in space and time that are very similar to other reports (Boman et al., 2018;Stoner and Appeldoorn, 2021).The derived percentages were then applied so that the number of particles released in the simulation follows the specified annual pattern (SI Table 1).

Unexploited and Exploited Spawning Scenarios
To account for the spatially heterogeneous patterns in exploitation and density-dependent processes (Allee effects), we considered two distinct scenarios.For the first set of simulations, the total number of larvae released from each polygon centroid is constant across all countries, with spawning varying over time to capture local seasonality.This set of simulations reflects a historical virgin or "unexploited" population, whereby the density of queen conch spawning biomass is constant across all areas.The actual densities of conch spawning biomass that can be supported will be dependent on habitat attributes of each area; however, we lacked a comprehensive map of habitat types across our wider study region, as well as information on the relationships between habitat types and their respective conch densities at carrying capacity, as they vary across the region.In the absence of such detailed information, we assumed that equal spawning biomass densities and consistent per-capita fecundity rate across the region are reasonable approximations for understanding overall connectivity patterns in a preexploitation (i.e., unexploited) scenario.
For the second set of simulations, densities of adult conch populations for each of the 42 jurisdictions of our study area (Horn et al., 2022, Figure 3) were used to scale the conch reproductive output.We applied an additional penalty for very low densities to account for reduced reproductive success at these levels (i.e., Allee effects).This simulation represents the expected connectivity under contemporary exploitation rates.To adjust the number of egg releases on a per-area basis, we used Equation 1, where N is the number of particles released in each year of the simulation per area a, S is a constant scalar to adjust the overall number of particles to computationally viable quantity for the simulation, D is the median of the mean adult conch densities for each jurisdiction c, and f(A) is a function representing the reduction in reproductive output due to Allee effects.
The few studies that have investigated simultaneously the relationship between adult density and reproductive rates found varying results, particularly because the density is dependent on the scale over which observations are recorded (Stoner and Ray-Culp, 2000;Goḿez-Campo et al., 2010;Stoner et al., 2012;Cala et al., 2013;Doerr and Hill, 2018;Appeldoorn, 2020).However, field studies suggest a cross-shelf minimum density of 50-100 adults/ha (Stoner and Ray-Culp, 2000;Stoner et al., 2012) or a spawning aggregation (i.e., intra-aggregation) density of 200-300 adults/ha (Delgado and Glazer, 2020) is necessary for successful encounters and mating.For our simulations, we assumed that jurisdictions with cross-shelf densities >100 conch/ha had full reproductive activity [f(A) from Equation 1 is set to 1], with reproductive output decreasing linearly for jurisdictions harboring between 50-100 conch/ha and assumed no reproduction for jurisdictions with densities <50 conch/ha (SI Figure 2).The number of particles released per polygon and jurisdiction was distributed across the annual cycle using the same derived percentages used in the unexploited, uniform spawning, simulation (SI Table 1).Applying this function to our contemporary estimates of adult queen conch density limited spawning to just 15 of the initial 42 jurisdictions (Figure 3 and SI Figure 4).

Vertical Migration Behavior
Hatching of veligers occurs from egg masses on the seafloor.Hatched veligers possess two lobes which are used to immediately swim up to the ocean surface and remain in these layers (Brownell, 1977).Because hatched veligers would reach the surface in a relatively brief amount of time, we did not simulate the process of hatching and swimming to the surface, instead simulating releases directly from surface waters where the early-stage larvae occur more frequently.
We did not simulate diel vertical migration since there is disagreement among studies as to whether significant patterns in diel distributions exist (Chaplin and Sandt, 1992;Stoner and Davis, 1997).However, the body of literature on conch early life history is largely in agreement that distributions shift as veligers become older, and thus we simulated this ontogenetic shift based on the comprehensive surveys from Stoner and Davis (1997).Stages of development were differentiated in the Stoner and Davis (1997) study as newly hatched (300-500 mm length), midstage (500-900 mm), and late-stage (>900 mm) larvae.To determine the timing of these phase changes in terms of time elapsed since hatching, we plotted reported ages and lengths from a suite of nine different studies (D'Asaro, 1965;Brownell, 1977;Boidron-Metairon, 1992;Davis et al., 1993;Davis, 1994;Davis et al., 1996;Davis, 2000;Brito-Manzano and Aldana Aranda 2004;Davis, 2005).The relationship between shell length and days post-hatching is linear and highly correlated (R 2 = 0.95).We used the relationship as defined by linear regression to translate shell length to estimated age at the transitional phases, which resulted in an estimated transition to mid-stage at 5.6 days post-hatching and transition to late stage at 15.5 days post-hatch.The initial vertical distributions (newly hatched veligers) are based on the averages of larval distributions with depth for newly hatched veligers from Tables 1, 2 in Stoner and Davis (1997).The vertical distributions of mid-and latestage veligers were extracted from Stoner and Davis (1997), where day and night vertical densities were averaged to determine an overall distribution.The data compilation showed that most early-stage larvae are found in the surface layer, with some distributed slightly deeper within the top 15 m, whereas the mid-stage and late-stage larvae are distributed slightly further down in the water column (from 5-30m), and to a lesser extent in the deeper layers (up to 100m, SI Figure 5).This summary is consistent with other reports of OVM behavior of veligers (Brownell, 1977;Barile et al., 1994;Davis and Stoner, 1994;Stoner et al., 1997;Stoner and Smith, 1998).

Settlement
Minimum pelagic duration is reported in three field studies to be 18 days (Brownell, 1977;Salley, 1986;Davis, 1994) and thus, we set the minimum duration to this value.Upper limits for pelagic larval duration range from 21 to 30 days (D'Asaro, 1965;metapoluation Brownell, 1977;Salley, 1986;Davis, 1994;Paris et al., 2008) with a mean of approximately 25 days; this value was used for the upper range of the allowable PLD in the simulation.

Post-Processing Analysis
Connectivity CMS saves individual larval positions and attributes through time, together with connectivity files with the spawning and settlement sites of successful settlers.The connectivity output is used to generate transition matrices of connectivity (i,j), where the rows (i) represent the spawning sites and the columns (j) represent the settlement sites (Bodmer and Cavalli-Sforza, 1968;Paris et al., 2007).By normalizing the connectivity matrix by columns, i.e., by the total settlement originating from each spawning site (i), we obtain the probability of settlers from a spawning site i to reach a settlement site j (Pij), and the matrix main diagonal (Pij, j=i) represents the probability that larvae settle within the same jurisdictional boundaries in which they were spawned (Paris-Limouzy, 2011).When the transition matrix is normalized by rows, or by the total settlement at each settlement site (j), the main diagonal of the matrix (Pji, i=j) represents self-recruitment, or the fraction of settlers spawned locally over all settlers in a reef (sensu Paris and Cowen, 2004;Botsford et al., 2009;Burgess et al., 2014;Drury et al., 2018).We also calculated the local retention, which is given by the number of settlers in their natal reef (or jurisdiction) over all larvae spawned in a reef (or jurisdiction, sensu Paris and Cowen, 2004;Botsford et al., 2009;Drury et al., 2018).Together these measures give an indication of how reliant any given jurisdiction is on larvae sourced from other jurisdictions.
Based on the transition matrix, we used graph theory to identify ecological corridors.Graph theory is the study of relationships between objects in a dynamic system, which can range from metapopulations to social or road networks (West 2001).A graph is composed of nodes and the connections between nodes are represented by edges.Thus, graphs are apt representations of the connections between populations created by larval dispersal (Cowen et al., 2006;Treml et al., 2008;Urban et al., 2009;Thomas et al., 2014;Saunders et al., 2016), where the populations are the "nodes'' and the connections between populations are the "edges."Several network metrics can be derived from graphs and here we calculate the betweenness centrality (BC), which measures the relative "centrality" of a node in relation to other nodes of the network.For each population (i), their BC is defined as the number of shortest paths (i.e., the shortest routes) between all population pairs which passes through the target population (i) over the number of shortest paths between all nodes in the graph.Thus, BC is a relevant measure of overall "connectedness," since the shortest paths include both direct and indirect -or multigenerational -connections (Treml et al., 2008;Holstein et al., 2014).Jurisdictions with high BC are generally important for maintaining genetic diversity -if reproductive output from jurisdictions with high centrality were to decline significantly, we would expect reduced genetic mixing over the region.Other metrics, such as in-degree and out-degree (number of connections towards and from a population, respectively), are also estimated.We used the Matlab-bgl package (Gleich, 2008) for all graph analyses.We also estimated a source-sink index (Holstein et al., 2014), which relates the number of settlers leaving or arriving at each jurisdiction's reefs (SS = settlers outsettlers in /settlers out + settlers in ).This index allows the characterization of relative sinks and sources within the region.
We calculated transition matrices using all the conch habitat polygons in the region (1359 polygons), and then grouped the polygons by jurisdiction of origin to calculate transition matrices and conduct other analyses.To create jurisdiction categories, we used boundaries from Exclusive Economic Zones (EEZ) from the Maritime Boundaries Geodatabase (Flanders Marine Institute, 2019).We defined population units within these geographical boundaries because: 1) both habitat and densities of queen conch were available by jurisdiction, and we are unable to provide robust input within jurisdictions and 2) these are the units of management throughout the queen conch distribution.To facilitate the description of results and discussion, the islands and jurisdictions of the Lesser Antilles (Windward and Leeward Islands, and the Leeward Antilles) are often described together, as color-coded and described in Figure 2.

Uncertainty Analyses
To account for the variable range of adult densities observed by fishery-independent surveys at each jurisdiction, we conducted an uncertainty analysis.From the above-mentioned jurisdictionspecific distributions of adult conch densities, we randomly drew 100 draws for each jurisdiction, and then translated the randomly drawn densities of adults to reproductive output as with previous simulations.We then estimated larval exchange between spawning and settlement jurisdictions by multiplying the probabilistic connectivity matrix of the unexploited release case by each reproductive output (egg production * connectivity, Kool et al., 2010;Kool et al., 2011).The results of the 100 transition matrices were aggregated for two extremes: minimum and maximum estimated connectivity.The matrices for each scenario were based on locating the minimum (or maximum) number of larvae exchanged between each spawning-settlement pair for the 100 different possible connectivity results.

Metapopulation
To identify clusters of highly isolated spawning and settlement sites of queen conch within the Wider Caribbean, we applied the graph theory method described in Jacobi et al. (2012).This method is a minimization algorithm that recursively seeks solutions that minimize connections between clusters, while maximizing connections within clusters (Jacobi et al. (2012), Garavelli et al., 2014).Briefly, it minimizes an objective function consisting of the connectivity between populations and a penalty parameter to avoid aggregating sites i.e., to avoid the optimal minimization solution of the connectivity matrix, which would be achieved considering only one cluster composed of all spawning and settlement sites.We used the Jacobi et al. (2012) minimization algorithm from the open-source R-package ConnMatTools (Kaplan et al., 2017).We did 10 different iterations of the algorithm using different penalty parameters (Jacobi et al., 2012).From these iterations, we selected the solution yielding the least amount of connectivity between clusters while resulting in the highest self-recruitment within the cluster.For this, we selected the solution where the connectivity between clusters normalized by the selfrecruitment resulted in the smallest values, a value termed "quality" (Jacobi et al., 2012).
Given the distribution of habitat (patches of habitat on coasts, islands, and banks) within the queen conch distribution area, we assumed that each one of these clusters, as characterized by the minimization algorithm and estimated larval dispersal, acts as a metapopulation (sensu Wells and Richmond, 1995).Thus, metapopulations were the results of the clustering analyses, populations were defined by jurisdictional divisions, and distinct banks within a population were considered subpopulations.

Estimation of Genetic Distance
From the probability matrices, we also derived the connectivity distance between each spawning-settlement pair for several generations.We did this by a forward projection of the matrix using a modified model of migration based on Bodmer and Cavalli-Sforza (1968) as described in (Kool et al. 2010;Kool et al. 2011) and Foster et al. (2012).In the unexploited scenario, we assumed that all habitats in our model have the same carrying capacity and reproductive output, and thus, the forward matrix with the allele frequencies expected at time t (Qt) is given by Q t = M + It, where M is the connectivity matrix, I is an identity matrix of the same dimension as M, and the bar denotes row-normalization (Kool et al., 2010).While the connectivity matrix reflects single dispersal events between populations during one generation, the forward projection of the matrix captures indirect connections over multiple generations.We then estimated the Nei's D A genetic distances (Nei's angular distance, as defined in Nei et al., 1983) between sites using the forward matrices, following Kool et al. (2010).The estimated Nei's D As are used for comparisons to the linearized G ST (G ST /(1-G ST )) of a Caribbean-wide genetic study (Truelove et al., 2017) by means of a Mantel-r test (Mantel, 1967;Legendre and Legendre, 2012) to test for isolation by distance using the R package adegent v2.1.1 (Jombart et al., 2008).

Connectivity Patterns
The probabilities of connectivity for both hydrodynamic fields (HYCOM and GLORYS, Figure 4) show that larval exchange tends to follow the main pattern of currents in the Caribbean, Gulf of Mexico, and North Atlantic (Figure 1).In the unexploited scenario, all network nodes (jurisdictions) were generally well connected, with multiple pathways of dispersal for both exporting and receiving larval subsidies -the median number of exporting connections (considering all jurisdictions and both hydrodynamic models) was 9, while the median number of incoming connections was 7. The connectivity from both hydrodynamic products also presented a high likelihood of self-settlement.However, the level of self-retention varied for some jurisdictions by model, e.g., Turks and Caicos Islands, which exhibited a higher likelihood of self-settlement when dispersed with the HYCOM than the GLORYS.Moreover, both networks showed that the only completely isolated node in the network is Bermuda, which did not receive nor provide larval subsidies to any other jurisdiction.
Some key differences between the two ocean circulation models were confined to regional scale connections.For instance, the HYCOM estimated a higher likelihood of exchange 1) between and within reefs in Cuba and The Bahamas, and 2) from Cuba to Mexico and the Greater Antilles, while the GLORYS simulations estimated a higher likelihood of export from Cuba to the Bahamas.Transport within the Windward and Leeward Islands also differed slightly between the two sets of simulations.For the HYCOM, larvae spawned on the Windward and Leeward Islands were estimated to have a higher likelihood of settlement in Puerto Rico and the Dominican Republic when compared to the high likelihood of local settlement observed in the GLORYS for this region.
When we introduced realistic reproductive output levels based on exploitation patterns (Figures 4B, D), connectivity was largely reduced among regions for both models, as expected, and differences in connectivity patterns between the two models became less apparent.Notably, when considering exploited densities (Figures 4, 6), there was a reduction in the likelihood of self-settlement.Indeed, even the highest probability of selfsettlement in the exploited scenario (The Bahamas, Jamaica, and Cuba) was smaller than observed in the historical scenario.Other jurisdictions in the network are likely to export most of their larval subsidies.By contrasting the connectivity networks for unexploited and exploited densities, there was a clear reduction of connections between the eastern and western Caribbean.It was also clear that the lack of spawning in Belize and Mexico, and the reduced spawning in Nicaragua, further break the transport to upstream sites, particularly the USA, The Bahamas, and northern Cuba.The regions that remained mostly interconnected were the Central American jurisdictions (Costa Rica, Honduras, Nicaragua), the westernmost Greater Antilles Islands (Jamaica, Cuba, Turks and Caicos Islands), and the offshore Serrana Bank.Larvae released in these regions were exported to other regional reefs (Belize, Mexico, Panama, Colombia, Haiti, Serrana, and Quitasueño Banks) for both models.
Generally, the pattern of export or import of larvae for the unexploited homogenous spawning scenario matched the experiments with both hydrodynamic fields (Figure 6A, "unexploited").Mexico, Costa Rica, British Virgin Islands, Anguilla, Saint Martin, and Saint Barthelemy presented a higher likelihood of acting as relative sinks, while the Windward Islands, Nicaragua, and Colombia's Mesoamerican offshore banks (Quitasueño, Serrana, and Serranilla) appeared as (relative) sources.Once reproductive output was scaled by adult conch densities (Figure 6A, "exploited"), the relative importance of regions as a source and a sink also changed.Saba, Saint Kitts and Nevis, Nicaragua, Serrana Bank, and Turks and Caicos Islands were sources on both models, while Puerto Rico, Honduras, and Cayman Islands became sinks.Two countries with the largest habitat areas, Cuba and Jamaica, appeared as neutral, which is expected since these countries also showed high self-recruitment.The probability of The Bahamas acting as a relative sink slightly increased in the exploited simulations.
The betweenness centrality (BC) for the unexploited scenario (Figure 6B, "unexploited") showed that despite differences in connectivity patterns between the simulations using different hydrodynamic products, the jurisdictions showing the highest BC were mostly consistent among the two simulations.These were the Dominican Republic, Puerto Rico, Jamaica, and Cuba.Following these were also Honduras, Colombia, and the British Virgin Islands.When considering exploited densities to scale the spawning (Figure 6B, "exploited"), Jamaica and Cuba emerged as central nodes for both hydrodynamic products.The Bahamas, Nicaragua, and Honduras appeared as secondary nodes.A notable change between the unexploited and exploited scenarios was Puerto Rico -when exploited densities were considered, the local reefs were unlikely to serve as an ecological corridor between Caribbean regions.
For the unexploited spawning scenario, local settlers comprised most of the settlers in several jurisdictions, including the Greater Antilles Islands, some of the Windward Islands (Figure 6C, "unexploited"), The Bahamas, Florida, and Bermuda, and results obtained with the two hydrodynamic products largely agreed.When considering exploited distributions (Figure 6C, "exploited"), local settlers comprised most of the settlers in Nicaragua, Jamaica, Cuba, and Turks and Caicos Islands, while largely important for The Bahamas and Saint Kitts and Nevis.Many of the jurisdictions that appeared to mostly rely on their own settlers in the unexploited experiment did not exhibit spawning (e.g., Haiti, Dominican Republic, Barbados, Trinidad and Tobago, Colombia, Belize, US, and Bermuda).

Uncertainty Analyses
Based on the minimum exchange scenario for each jurisdiction (Figures 5A, C), connectivity was further reduced from the exploited scenario, with more countries not producing viable settling larvae (or not spawning), including Honduras, Turks and Caicos Islands, Serrana Bank, St. Kitts and Nevis, and Sint Eustatius.However, the reduction in spawning did not further alter the probability of connectivity with other jurisdictions which were spawning, either for the HYCOM or GLORYS simulations.However, there were differences in the strength of connections and the number of larvae settling on their natal reefs.There was a higher likelihood of self-settlement in Jamaica, Cuba, and The Bahamas for both models.The connectivity between the eastern and western Caribbean regions here was further fragmented.
On the other hand, for the scenario considering the maximum connectivity between two countries (Figures 5B,  D), spawning occurred in nine additional jurisdictions from the sixteen in the exploited scenario, including the British and US Virgin Islands, Guadeloupe, Colombia, Serranilla Bank, Belize, Mexico, and The Bahamas (Western Great Bahama Bank and Central Great Bahama Bank).As expected, this scenario presented higher connectivity than the exploited scenario (based on median values of densities), although a break in larval transport between the western and eastern Caribbean still emerged, since the centrally located Dominican Republic and Haiti only served as settlement sites and did not support reproductively active populations in this scenario.Several countries also presented a higher likelihood of local settlement in this case, particularly Jamaica, Cuba, The Bahamas, and Colombia.Patterns of transport and local settlement differed between the two hydrodynamic models for the Windward and Leeward Islands, but generally there was a higher likelihood of local settlement for the jurisdictions still reproductively active in this case.The probability of connectivity was still similar for the countries where conch were spawning in the two scenarios.However, in this case, there was much more exchange between the Western Caribbean jurisdictions and a higher likelihood of export of larvae downstream of the local currents, particularly to countries located downstream, such as Mexico and Florida.

Metapopulation
Regional patterns of connectivity between distinct Caribbean and North Atlantic regions showed changes between hydrodynamic models and scenarios (Figure 7 and SI Table 3).The major differences in the unexploited scenario between models were estimated in Cuba, south Florida, the Leeward Antilles, southern Mesoamerica (Costa Rica and Panama), and the southern Leeward Islands for the unexploited scenario and both hydrodynamic models, with six metapopulations.For the HYCOM simulations there were: 1) the Lesser Antilles and Puerto Rico; 2) Venezuela, the Leeward Antilles, northwestern Caribbean, Florida; 3) Haiti, Dominican Republic, The Bahamas, Turks and Caicos Islands; 4) Cuba and Cayman Islands; and 6) Bermuda.For the GLORYS: 1) Western Caribbean, Jamaica, and Cuba; 2) Florida; 3) Puerto Rico and Northern Lesser Antilles; 4) Southern Lesser Antilles; 5) The Bahamas, Turks and Caicos Islands, and Venezuela; and 6) Bermuda.
When spawning was reduced according to the relative exploited densities of conch, connections were decreased by the lack of spawning in exploited jurisdictions, as particularly evident in Florida, Mexico, Haiti, and the Dominican Republic.However, the entire region was still connected by larval dispersal through connections within the central regions of the network.When exploited densities were considered, there was an increase of metapopulations for the HYCOM (from 6 to 7 metapopulations, SI Table 3) and a decrease for the GLORYS (from 6 to 4 metapopulations).Remarkably, for both models, the jurisdictions harboring the largest conch populations (Cuba, Jamaica, and The Bahamas) became metapopulations by themselves (Jamaica) or only with their nearest neighbors (Cuba and Cayman Islands, The Bahamas, and Turks and Caicos Islands).Both models also diverged in the number of connections exhibited by each node.
The metapopulations obtained with the uncertainty analyses (SI Figure 6) suggested that the metapopulation organization of queen conch is highly dependent on the spatial patterns of reproductive output.The metapopulations obtained with the minimum connectivity (SI Figures 6A, C) were highly distinct for the two hydrodynamic models.While six connected metapopulations were identified for the HYCOM, only two isolated metapopulations emerged for the GLORYS -one encompassing the western Caribbean, and another for the easternmost region (Windward and Leeward Islands and Puerto Rico).For the maximum exchange, the number of metapopulations remained the same for the HYCOM, although with changes to its spatial structure.The number of metapopulations greatly increased for the GLORYS, which showed eight distinct clusters of mostly self-contained habitat.Horn et al., 2022).Minimum (maximum) potential connectivity is obtained from an uncertainty analysis, where we multiply the probabilistic connectivity matrix of the unexploited release case by 100 spawning distributions drawn from the density curves of adult conch per jurisdiction and use the smallest (biggest) number of larvae exchanged between each spawning-settlement pair to generate the probabilistic connectivity.Edges between nodes represent the strength of connections between sites, whereby the direction of larval flow is given by the clockwise direction -concave edges represent export while convex represents import.Node sizes represent the probability of self-settlement.
Metapopulations were more similar among the models, i.e., Panama, The Bahamas, and Jamaica were individual metapopulations for both models.Similarly, output from both models resulted in a cluster formed by jurisdictions not exhibiting reproductive activity and that showed a low probability of receiving larval subsidies from other metapopulations.

Ecological and Genetic Connectivity
The origin of settlers at each discrete habitat unit was grouped by the genetic divisions characterized by Truelove et al. (2017, Figure 8).The unexploited spawning scenarios (Figures 8A, B) corresponded well to the patterns observed by Truelove et al. (2017) given that larvae within each region identified by the authors were most likely locally originated.The exception was the high probability of larval exchange between The Bahamas and Turks and Caicos Islands and the Greater Antilles, which occurred for both hydrodynamic models.In the exploited reproductive output scenario, we considered that six of the twelve jurisdictions sampled by Truelove et al. (2017) were not reproductively active.Due to the lack of spawning, it was expected that not all connectivity patterns could be reproduced.Indeed, in this case, the high self-settlement observed for Mexico, Belize, and Florida (Figures 8C, D) was absent due to the lack of reproductive activity.The genetic connectivity analyses subsequently focused only on the results of the unexploited simulations since the results of the exploited density scenario were insignificant due to the reduced number of data points.
The results of the Mantel-r test between G ST values from Truelove et al. ( 2017) and the Nei's D A genetic distance derived from our connectivity matrices for 10, 50, and 100 generations (SI Table 2) showed that the Mantel-r increased with the number of generations considered.By considering the averages of the Nei's D A genetic distances calculated for each individual hydrodynamic product, the Mantel-r increased for 50 and 100 generations, while the p decreased.This approach captured the uncertainties of the hydrodynamic field estimated by each individual model and improved our estimations (Figure 9 and SI Figure 7).Results of both the Mantel-r tests (SI Table 2), the scatter plot between the estimated Nei's D A genetic distances and linearized G STs from Truelove et al., (2017, Figure 9), and plots of distances between sites (SI Figure 7) suggested that queen conch populations exhibit a distance-by-isolation pattern.Our model explained part of the variability and the consideration of conch specific biological traits improved the results obtained solely by considering oceanographic distances between sites.

DISCUSSION Connectivity Patterns
The connectivity network (Figure 4) of queen conch for both hydrodynamic models used generally follows the main surface currents in the region -Caribbean Current, Yucatań Current, Loop Current, and Florida Current (Figure 1) -flowing in a clockwise direction from the Lesser Antilles to the Florida Straits.In the northern part of the Greater Antilles, the Antilles Current flowing northwestward favors connections within the Windward and Leeward Islands with the Greater Antilles and Bahamas.Despite the remarkable footprint of this current system, the circulation in the Inter-American Sea is highly variable.The intensity and extension of the currents change at seasonal and inter-annual scales and numerous sub-mesoscale and mesoscale features facilitate recirculation, settlement of larvae close to their spawning sites, and larval exchange with neighboring jurisdictions (Paris and Cowen, 2004;Paris et al., 2005;Kough et al., 2013;Kough et al., 2016;Lopera et al., 2020).Indeed, numerous recurrent gyres are known to occur in the region, notably the semi-permanent cyclonic Panama-Colombia Gyre (Lopera et al., 2020), the transient Honduras and Cozumel Gyres (Carrillo et al., 2015), and the anti-cyclonic gyre around Cayman Islands (Molinari, 1980).Thus, differences in connectivity between the hydrodynamic products reflect the variability of these surface and subsurface flow patterns, particularly around the Lesser and Greater Antilles, which affect the connectivity patterns between The Bahamas, Cuba, and the Windward and Leeward Islands.However, the most apparent differences in connectivity are between our exploitation scenarios, and not between our hydrodynamic products.
The most apparent differences between the unexploited and exploited connectivity emerged from the fact that most jurisdictions were estimated to currently hold conch densities well below the critical threshold for reproduction (Figure 3, Horn et al., 2022), and thus were here considered to be reproductively non-functional (SI Figure 4).These jurisdictions, therefore, could only act as sinks for larvae, not providing larval subsidies for other locations or themselves.The connectivity patterns emerging from the exploited scenarios, considering any range of densities (maximum,minimum,or median Figures 4,5), are thus drastically different than expected from uniform spawning in the unexploited scenario.For example, due to their up-current position, the Leeward and Windward Islands (i.e., Lesser Antilles) were estimated to be historically important for contributing larval input to other jurisdictions, particularly the Leeward Antilles, Venezuela, and the Greater Antilles.However, due to overfishing, most of these jurisdictions are no longer likely to provide larval subsidies in the FIGURE 7 | Connectivity network analyses for queen conch metapopulations identified with a minimization algorithm (Jacobi et al., 2012;Garavelli et al., 2014), for two different reproductive output scenarios: 1) unexploited, i.e., uniform spawning for all habitat polygons (A, C left column); and 2) exploited, or spawning scaled by median adult densities reported per jurisdiction (B, D right column).Metapopulations are identified by the habitat polygon colors (gray habitat represents nonreproductive populations, i.e., where no spawning occurs).Larvae are dispersed with HYCOM velocities (A, B, top row), and GLORYS12v1 reanalysis (C, D, bottom row).Edges between nodes represent the strength of connections between regions, whereby the direction of larval flow is given by the clockwise direction -concave edges represent export while convex represents import.Node sizes represent the probability of self-settlement -their location is given by the mean latitude and longitude of all habitat polygons in that specific metapopulation.Colors vary between scenarios, since the number and configuration of metapopulations also change.
exploited scenarios and are relatively isolated from outside replenishment.
Other patterns in comparing the unexploited and exploited simulations were more subtle, but locally significant.For example, historically the Turks and Caicos Islands were estimated to have received subsidies from the Dominican Republic and Haiti.However, due to overfishing in these source jurisdictions, the exploited scenario suggests that the Turks and Caicos Islands are now dependent on reduced local production, which is less likely to be locally retained, with a large percentage of larvae being exported to The Bahamas.Likewise, the unexploited simulation suggests that Florida was dependent on relatively high local retention, with other significant input coming from Mesoamerica (particularly Nicaragua, Belize, and Mexico).These jurisdictions now have estimated densities below the critical reproductive threshold; thus, Florida currently has no significant sources of larvae, except a small likelihood of receiving larvae from Cuba.This could explain why, despite a moratorium on fishing for several decades, Florida's conch population has not recovered (Delgado and Glazer, 2020).
Heterogeneous fishing pressure and localized depletion also significantly influenced the role of each jurisdiction in the conch network (Figure 6), as exemplified by the changes in the betweenness centrality of the Dominican Republic and Puerto Rico (Figure 6B).These jurisdictions presented some of the highest centrality measures in the unexploited scenario yet no longer function as important connectors in the exploited scenario.An almost complete break in the network was apparent based on the exploited scenario, with the Dominican Republic receiving limited larvae from Cuba, Turks and Caicos Islands, and from a deep mesophotic reef off the west coast of Puerto Rico.The channel between Puerto Rico and the Dominican Republic, the Mona Passage, has been previously shown to act as a transport barrier or filter between these two jurisdictions (Baums et al., 2006;Beltrań et al., 2017).However, the flow in the channel is highly variable (Metcalf et al., 1977), allowing for a small likelihood of periodic connections between populations across the passage (Baums et al., 2006;Beltrań et al., 2017;Prakash et al., 2018, Willis et al., 2022).Since our putative mesophotic spawning site was located at the Abrir La Sierra reef (Garcia-Sais et al., 2012), on the eastern edge of the Mona Passage, it is expected that probabilities of cross-channel connections might be observed in the explored scenario.
While the Dominican Republic and Puerto Rico nodes have been removed from the chain of larval supply, Jamaica and Cuba remained as important connectors in the western portion of the range, and some of the offshore banks in Colombia remained as functional connectors.Although connections have been lost in other locations due to localized overfishing, there were still connection points, albeit reduced, that would allow some exchange of larvae and maintenance of some genetic diversity over time.

Ecological and Genetic Connectivity
Localized patterns of conch overfishing also affected the conch network in the Wider Caribbean Region and could influence both conch demographics and genetic population structure.As shown by the clustering analyses of sub-populations under unexploited and exploited reproductive output (Figure 7 and SI Figure 6), changing the reproductive output dramatically In the unexploited scenario, individual settlement habitat units (i.e., polygons) in Mexico and Belize showed a high level of self-recruitment (Figure 6C), which corroborated the genetic patterns observed by Truelove et al. (2017).However, when considering larval exchange probabilities at a jurisdiction level, this pattern was not clear and the populations were not isolated, as also demonstrated by the metapopulations emerging from the metapopulation analyses (Figure 7 and SI Figure 6).Furthermore, these patterns disappear when considering exploited contemporary density estimates of adult conch, highlighting a disconnect between genetic and ecological (or demographic) connectivity of this population.This disconnect is not unexpected since genetic and ecological connectivity explore processes spanning different spatiotemporal scales.The ecological connectivity covers the shorter scales of early planktonic stages, examining the exchange of individuals between populations (Kough et al., 2017).Typical studies of genetic connectivity unveil patterns occurring through multiple generations or even evolutionary timescales and characterize the gene flow between these populations (Thomas et al., 2014;Brown et al., 2017;Carr et al., 2017).Whereas simulations of larval dispersal can inform the evolutionary connectivity of organisms, characterizing, for instance, transport barriers and isolated populations (Cowen et al., 2006;Baums et al., 2006;Truelove et al., 2017), dispersal models are particularly useful to help elucidate patterns of ecological connectivity (Jones et al., 2009).Indeed, our results reflect the different spatial-temporal scales unveiled by the genetic and ecological connectivity.By comparing similarities for multiple generations obtained from the unexploited scenarios to genetic divergence (linearized G ST ), we found these to be significantly correlated.However, once we considered localized depletion patterns, these correlations disappeared, suggesting that the genetic results reported by Truelove et al. ( 2017) could have captured patterns of connectivity that no longer exist under the (conservative) depletion scenarios considered in our study.

Management Implications
The objective of this study was to explore how localized harvest patterns affect the probabilistic ecological connectivity of the queen conch throughout its range.The lack of detailed demographic and habitat information (e.g., population distribution, age structure, adult densities, fecundity, carrying capacities) precludes a direct estimation of the strength of larval recruitment.Regional cooperative management efforts have advocated for a more consistent and uniform collection of these data across the WCR (Prada et al., 2017).These data would allow us to calculate the persistence and replenishment rates of individual populations (Leis, 2006;Burgess et al., 2014;Lowerre-Barbieri et al., 2017;Hamilton et al., 2021).Moreover, the interpretation of ecological and genetic connectivity has important spatial management implications.While small and sporadic exchanges of a few individuals might lead to panmixia of distant sub-populations, a population requires a steady and strong source of larval subsidies for its support (Wright, 1931;Botsford et al., 2003;Lowe and Allendorf, 2010;Paris-Limouzy, 2011;Hawkins et al., 2017).Therefore, ecological connectivity projections provide vital information for management by delineating areas that are self-reliant on their recruits to support their populations and in identifying relative sources and sinks (Botsford et al., 2003;Botsford et al., 2009;Hawkins et al., 2016).
Here, we characterized distinct metapopulation delineations for queen conch for the scenarios considered in our study (distinct hydrodynamic models and assumptions of adult densities) throughout its range.The existence of these self-sustained metapopulations implies the need for effective management measures within these regions.Here we showed that configuration of such metapopulations is strongly affected by localized depletion patterns, more so than by oceanographic patterns, given the changes in the metapopulation configurations in the homogenous and depleted scenarios (Figure 7).As these metapopulations cover multiple jurisdictions, their management would thus require international cooperation (Kough et al., 2017).
Connectivity processes have been explored in many studies using larval dispersal simulations in the last decades; however, studies still consider uniform spawning populations, despite the noted importance of the effects of variable spatio-temporal fecundity on connectivity patterns (Hughes and Tanner, 2000;Karnauskas et al., 2011;Lowerre-Barbieri et al., 2017;Castorani et al., 2017;Johnston et al., 2018;Hamilton et al., 2021).We found that the inclusion of realistic present-day patterns in reproductive output changes connectivity patterns and the structuring of queen conch sub-populations and metapopulations across their range, which must be considered in future studies.
Overharvest not only can change the patterns of connectivity, as addressed in this study, but also lower the genetic diversity of fisheries populations (Pinsky and Palumbi, 2014).A reduction in diversity due to conch overharvest has been suggested by Zamora-Bustillos et al. (2011) investigating genetic diversity and structure of two populations along the Yucatań peninsula in Mexico, and by Kitson-Walters et al. (2018) studying the genetic differentiation within the Jamaican EEZ.Truelove et al.
(2017) also observed a reduction in allele frequencies in Mexico and Anguilla, however, it is unclear whether these changes are due to isolation by distance, regional overharvest, local overharvest, or their combined effect.The results of our theoretical "virgin stock" simulations, however, suggested that there are no complete barriers to conch larval transport in our study area and thus, exploited patterns of isolation and inbreeding could be related to anthropogenic actions, such as habitat destruction and overharvest.This is a topic which warrants further investigation.
The simulations relied on several assumptions regarding the spatial structure and demographics of conch populations in the Caribbean.As illustrated by our results (Figures 4-7), estimated ecological connectivity patterns are sensitive to assumptions regarding both the distributions of queen conch populations and their reproductive output.Thus, we highlight the need for updated demographic and habitat information for queen conch, strongly corroborating needs identified in previous studies (Prada et al., 2017;Kough et al., 2017;Truelove et al., 2017;Kitson-Walters et al., 2018;Boman et al., 2021).To refine larval connectivity estimates and fully understand conch demographics at both broad and smaller spatial scales within jurisdictionsthe following are needed: 1) representative surveys of conch population densities across a range of habitats, along with a precise understanding of the areas over which these densities can be extrapolated (Prada et al., 2017;Kough et al., 2017); 2) densities of individuals categorized as juveniles and adults and/ or detailed stock size structure reported; 3) detailed maps and/or areal estimates of conch habitat (i.e., areas where the species currently or has historically occurred) in each jurisdiction; and 4) surveys of conch populations located in depths below the typical range of fishing (i.e., >30m) to assess the contribution of mesophotic populations that might serve as deep refugia for shallow populations and which remain largely excluded from conch surveys (Boman et al., 2021).
Overall, our results support that estimates of contemporary demographic rates, and thus of fecundity, need to be considered for ecological connectivity studies (Lowe and Allendorf, 2010;Castorani et al., 2017;Johnston et al., 2018).Moreover, our results strongly corroborate advice for international management and conservation of queen conch, as well as research needs identified in previous conch studies (Prada et al., 2017;Kough et al., 2017;Stoner et al., 2019;Delgado and Glazer 2020;Boman et al., 2021).However, given the high likelihood of local settlement and distribution of distinct metapopulations, replenishment processes through larval dispersal are likely to occur only within the metapopulation spatial scales.Our simulations showed that local sub-population success would depend on effective local management, but also on the management of specific upstream source populations.The combination of local self-seeding and distant larval transport are both required for success throughout the species' range.Regional management efforts are working to implement many of the management plans and measures detailed by Prada and collaborators (2017).Local measures of conservation and management are thus critical to protect individual reproductively active sub-populations and nursery habitat across multiple jurisdictions.

FIGURE 1 |
FIGURE 1 | Average surface velocity fields (m s -1 ) and eddy kinetic energy (EKE, m 2 s -2 ) for HYCOM reanalysis and analysis (A) and GLORYS12v1 reanalysis (B).Figures are based on surface velocity fields from 2013 to 2017 (March to December), the same time period used for the simulations of dispersal.A span of 5 years was chosen to capture the interannual variability in the circulation and derived connectivity patterns.Arrows represent average surface velocity and background color the EKE.Labels in the top panel identify major oceanographic features in the region (CC, Caribbean Current; PCG, Panama-Colombia Gyre; YC, Yucatań Current; LC, Loop Current; FC, Florida Current; GS, Gulf Stream; AC, Antilles Current).

FIGURE 2 |
FIGURE 2 | Map of spawning and settlement habitat (polygons), color-coded by jurisdiction.Conch habitat is represented of 1359 8×8 km polygons (Holstein et al., 2014) based on shallow (depth less than 20 m) coral reefs locations from the Millennium Coral Mapping Project (Spalding et al., 2001; IMaRS-USF 2005; IMaRS-USF and IRD 2005; Andrefouët, 2008; UNEP-WCMC et al., 2021).Individual habitats are buffered to facilitate visualization and are not represented to scale.Due to close proximity, some jurisdictions are referred by their regional location, as follows: 1) the Greater Antilles is composed by the Cayman Islands, Cuba, Dominican Republic, Haiti, Jamaica, Puerto Rico; 2) the Leeward Islands are Anguilla, Antigua and Barbuda, British Virgin Islands, United States Virgin Islands, Guadeloupe, Montserrat, Saba, Saint Bartheĺemy, Collectivity of Saint Martin, Sint-Eustatius, Saint Kitts and Nevis; 3) the Windward Islands are composed by Barbados, Dominica, Grenada, Martinique, Saint Lucia, Saint Vincent and the Grenadines, Trinidad and Tobago.

FIGURE 4 |
FIGURE 4 | Connectivity networks for queen conch for two different reproductive output scenarios: 1) unexploited, i.e. spawning uniform for each habitat polygon (A, C left column); and 2) exploited, with spawning scaled by median adult densities reported per jurisdiction (B, D right column).Larvae were dispersed with HYCOM (A, B, top row), and GLORYS (C, D, bottom row) velocities.Edges between nodes represent the strength of connections between sites, whereby the direction of larval flow is given by the clockwise direction -concave edges represent export while convex represents import.Node sizes represent the probability of self-settlement.Jurisdictions with no spawning in the exploited density scenario are marked with a triangle.The results for each model were based on 90 releases in the period of five years (2013-2017).

FIGURE 5 |
FIGURE 5 | Connectivity network for queen conch considering uncertainty of contemporary (exploited) conch densities, for larvae dispersed with HYCOM (A, B, top row) and GLORYS12v1 (C, D, bottom row) velocities.Panels exhibit the minimum (A, C, left column) and maximum (B, D, right column) potential connectivity based on possible ranges of adult conch densities (fromHorn et al., 2022).Minimum (maximum) potential connectivity is obtained from an uncertainty analysis, where we multiply the probabilistic connectivity matrix of the unexploited release case by 100 spawning distributions drawn from the density curves of adult conch per jurisdiction and use the smallest (biggest) number of larvae exchanged between each spawning-settlement pair to generate the probabilistic connectivity.Edges between nodes represent the strength of connections between sites, whereby the direction of larval flow is given by the clockwise direction -concave edges represent export while convex represents import.Node sizes represent the probability of self-settlement.

FIGURE 6
FIGURE 6 | (A) Relative Source-Sink (SS) Index for each jurisdiction in the queen conch network.Positive values of SS Index represent relative source areas while negative represent relative sink areas.(B) Relative betweenness centrality (BC) for each jurisdiction in the queen conch network.BC measures the fraction of shortest paths passing through a node (jurisdiction) -higher values represent most central nodes, or ecological corridors.(C) Probability of self-recruitment for each jurisdiction in the queen conch network, defined as the proportion of settlers spawned locally in each jurisdiction (local settlers over all settlers in a given jurisdiction).Simulations considered two different reproductive output scenarios: 1) unexploited, i.e., uniform spawning for all habitat polygons; and 2) exploited, with spawning scaled by median adult densities reported per jurisdiction.Larvae dispersed with the HYCOM and the GLORYS12v1 velocities.When considering exploited density of adults, areas with no spawning are colored gray on a).Each node represents all polygons from individual jurisdictions, organized in a clockwise direction from Haiti towards Bermuda.Each model's results are based on 90 releases in the period of five years (2013-2017).

FIGURE 8 |
FIGURE 8 | Origin of larvae per settlement polygon, for two different reproductive output scenarios: 1) unexploited, i.e., spawning uniform for all habitat polygons (A, B); and 2) contemporary (exploited), with spawning scaled by median adult densities reported per jurisdiction (C, D).Larvae dispersed with HYCOM (A, C), and GLORYS12v1 (B, D) velocities.Origin areas are relative to genetic divisions identified by Truelove et al., (2017).

FIGURE 9 |
FIGURE 9 | Scatter plot of Nei's D A genetic distance (from simulations) and pairwise genetic differentiation (linearized G ST , based on G ST from Truelove et al., 2017) for Aliger gigas.The simulated Nei's D A genetic distance is an average of the genetic distances obtained with the HYCOM and GLORYS12v1 velocities for (A) 50 generations and (B) 100 generations.The trend line is shown in blue, and dots are color-coded by the sampling jurisdiction (or the source jurisdiction for the larval dispersal model).