Demographic connectivity links local populations via the dispersal of individuals in their larval, juvenile, or adult phases (Cowen and Sponaugle, 2009; Kendrick et al., 2017). In marine fish, the spatial scale of larval dispersal is generally orders of magnitude larger than in the juvenile or adult phases. Thus the transport of larvae via ocean currents is the main driver of the biological connectivity of many marine fishes (Green et al., 2015). Most studies of genetic connectivity in fish focus on sampling adult individuals, typically representing multiple cohorts and generations of larval and adult dispersal and after the effects of survival during settlement and growth, integrating all these processes into a snapshot of genetic structure and connectivity (Hedgecock et al., 2007; Gagnaire, 2020). Population genetic studies of larvae in transit are relatively rare. Still, they can shed light on an important yet complex stage in the life history of marine fish that eventually influences the population’s genetic process (Ottmann et al., 2016; Shima and Swearer, 2016; Riquet et al., 2017). Three distinct, non-mutually exclusive processes are in play to define the role of a site within a fish metapopulation, including the export of individuals towards other sites, the import of individuals from other regions, and the amount of self-recruitment (individuals that remain at their natal site at the end of their larval stage). The balance of these processes defines a site as a net source (export of individuals exceeds imports) or net sink (import of individuals exceeds exports), while the levels of exchange of individuals with other sites define if a population is considered open (high exchange) or closed (little or zero exchange) (Cowen and Sponaugle, 2009).

The availability of genetic information from larval pools is logistically challenging. The composition and fate of the larval pool itself depend on different physical and biological factors in the environment (Pineda et al., 2007; Otwoma et al., 2018). The environmental unpredictability is thought to cause a large proportion of larvae to die or eventually fail to settle and survive, with only a small subset of larvae surviving and eventually reproducing as adults (Shima and Swearer, 2016). This phenomenon is dubbed sweepstakes reproductive success and effectively reduces the effective population size and enhances genetic drift and the differences in allele frequencies among subpopulations (Ottmann et al., 2016). However, some studies have contradicted the sweepstakes hypothesis by showing genetic homogeneity among multiple larval events (Domingues et al., 2011).

The mixing process in the dispersal stage is essential for larval pools to maintain high genetic diversity (Yearsley et al., 2013). However, it has been reported that in some species, larvae can avoid the mixing process and display a collective dispersal. This implies that the movement of larvae is not entirely independent of one another (Yearsley et al., 2013). In the scenario of collective dispersal, larvae of the same spawning event tend to travel together, either driven by larval behavior or physical oceanographic processes. Collective dispersal allows relatives to stick together during the dispersal phase, increasing measures of endogamy in larval pools (Riquet et al., 2017).

The Twospot flounder (Bothus robinsi) is a small fish (max. length 30-40 cm) that inhabits the continental and insular shelf between depths of 0-100 m. In the Gulf of Mexico (GOM), B. robinsi is found from Florida to the Quintana Roo continental shelf (Munroe, 2003; Evseenko, 2008; McEachran, 2009; Blanqueto-Manzanero and Vega-Cendejas, 2018) ¹ .

Due to its small size, this flounder does not have a targeted commercial fishery. However, it has ecological significance because it contributes to the nutrient flux, acting as predator and prey (in all their developmental stages). They feed on small fishes, crustaceans such as crabs, shrimps, amphipods, polychaetes, among others organism. Bothus robinsi eggs and larvae are consumed by jellyfish, ctenophores, worms, shrimps, and fish during their time in the plankton. Young flounders are preyed by crabs, shrimps, and many predatory fishes (cod, sculpin, and moray eels, among other fishes) (Wakida-Kusunoki et al., 2013) ^{2 , 3} . Furthermore, B. robinsi larvae have the potential to be a connectivity bioindicator because of their high abundance and distribution in neritic and oceanic waters and their occurrence in all seasons (Houde et al., 1979; Sanvicente-Añorve et al., 2000; Flores-Coto et al., 2009).

Its principal dispersal strategy depends on the presence of pelagic eggs and larvae. Bothus robinsi has a high dispersal potential due to its capability of spawning throughout the year and its long pelagic larval duration (PLD) (Houde et al., 1979; Evseenko, 2008). Like other flounders, species from the Bothus genus can extend their larval developmental stage. A study of B. podas indicates that its larvae could settle after three months, depending on the environmental characteristics (Morato et al., 2007). For B. robinsi, there is a record suggesting the possibility for PLD as long as 5-12 months (Evseenko, 2008). Flatfishes from the Order Pleuronectiformes are characterized by long PLD, for example, 70 days for the Bay whiff (Citharichthys spilopterus; Gunther, 1862) (Joyeux et al., 1995), and 4-7 months for the Atlantic halibut (Hippoglossus hippoglossus; Linnaeus, 1758) (Sadorus et al., 2021).

Bothus robinsi dispersal potential does not end after finishing their larval stage. Their juvenile stage can live sometime in the pelagic environment until settlement because they use the reserve of nutrients in their massive liver to avoid starvation. In addition, they keep a transparent and slightly visible coloration that allows them to blend in the pelagic environment and avoid predators, increasing their probability of survival and settling (Evseenko, 2008). Bothus robinsi dispersal capabilities translate into a high potential for population connectivity in the GOM.

The GOM is a semi-closed sea with a 1.6 million km² area. This gulf is a dynamic system due to the presence of heterogeneous currents that fragment the pelagic environment into provinces characterized by some degree of isolation ( Figure 1 ) (Miron et al., 2017). The provinces were determined using lagrangian analysis based on the dynamical system that is governed by the motion of fluid particles, which are described by satellite-tracked drifters on the ocean surface (Miron et al., 2017). Starting the analysis, the tracers spend more time moving and eventually temporarily accumulating in regions within them than dispersing across their boundaries ( Figure 1 ). These retention effects, which constrain connectivity between distant places in the GOM and keep planktonic organisms inside the provinces, are the result of oceanic dynamic systems such as cyclonic and anticyclonic eddies, and circulation ( Figures 1 , 2 ) (Díaz-Flores et al., 2017; Miron et al., 2017).

The principal oceanographic feature in the GOM is the Loop Current and its related eddies. This current enters into the GOM through the Yucatan Channel, turns anticyclonically to the east, and leaves the GOM through the Straits of Florida ( Figure 2 ) (Díaz-Flores et al., 2017; Lara-Hernández et al., 2019). Anticyclonic eddies can shred any month of the year from the main Loop Current as a result of multiple seasonal and stochastic processes. These eddies drift toward the west or northwest ( Figure 2 ) until their arrival on the continental shelf between the north of Veracruz and the south of Texas, after 2-4 months, and dissipate (Muller-Karger et al., 2015; Díaz-Flores et al., 2017; Romero-Centeno et al., 2021). These large anticyclonic eddies (200-300 km) interact with the continental shelf, enhancing coastal upwelling throughout the region. The interaction of cyclonic and anticyclonic eddies allows the transport of water and particles among different provinces (Muller-Karger et al., 2015; Lara-Hernández et al., 2019; Romero-Centeno et al., 2021).

Although there is an overall anticyclonic circulation in the GOM, there are zones with cyclonic circulation, including the continental shelf between Texas and Louisiana, the shelf off western Florida and the cyclonic eddy in the Bay of Campeche ( Figure 2 ). Over the western continental platform (i.e.,Tamaulipas and Veracruz), seasonal winds drive the direction of the currents northward during May-August and southward or south-eastward during September-March ( Figure 2 ) (Lara-Hernández et al., 2019; Romero-Centeno et al., 2021).

Several studies in the GOM have modeled Langranian trajectories of particles (virtual larvae of propagules) driven by ocean currents using the Hybrid Coordinate Ocean Model (HYCOM) (Johnson et al., 2013; Sanvicente-Añorve et al., 2014; Lara-Hernández et al., 2019). Typically, passive particles are tracked during a period based on the species’ PLD to estimate the number of particles that settle in distinct regions and the patterns of potential connectivity (Johnson et al., 2013; Sanvicente-Añorve et al., 2014; Lara-Hernández et al., 2019). According, to simulations of larval transport (Johnson et al., 2013; Salas-Monreal et al., 2018; Lara-Hernández et al., 2019), the PLD influences the probability of larvae escaping their natal oceanographic province. The Campeche province is inside the Bay of Campeche, a southern bight surrounded on three sides by the Mexican states of Campeche, Tabasco, and Veracruz. This province has a tendency to retain plankton as a result of the seasonal neritic dynamics, the cyclonic eddy, and the influence of an anticyclonic eddy located to the north over deep waters ( Figure 2 ) (Sanvicente-Añorve et al., 2014; Villegas-Sánchez et al., 2014; Díaz-Flores et al., 2017; Lara-Hernández et al., 2019).

Studies modeling larval transport in the GOM suggest that fish with short PLD (16 days) have a high degree of connectivity among nearby subpopulations separated by ~56 km inside the Bay of Campeche. However, subpopulations located in the shelf of Veracruz and Campeche that are separated by ~460 Km exhibit lower connectivity (Sanvicente-Añorve et al., 2014; Lara-Hernández et al., 2019). Species with longer PLD (31 days) display potential connectivity between subpopulations from distinct provinces, including Campeche and Southwest Florida Shelf (SWFS), separated by 1280 Km (Johnson et al., 2013; Lara-Hernández et al., 2019). Organisms with a PLD 4-7 months like the lobster (Panulirus argus) show connectivity between distant subpopulations from different provinces separated by 1,680 km (Northern Veracruz and Southern Florida) (Lara-Hernández et al., 2019).

In this study, we contrasted allele frequencies between larval pools from multiple locations around the GOM and adult fishes of B. robinsi from the northern GOM. We also evaluated the genetic diversity in the larval pool of B. robinsi and tested if the larvae tended to mix or if they displayed evidence of collective dispersal. Finally, we evaluated the genetic connectivity using B. robinsi larvae and adults across the GOM.

A total of 105 larvae were molecularly identified as B. robinsi ( Table 1 and Supplementary Table 2 ) and 21 larvae were molecularly identified as B. ocellatus. Overall, most identified samples were at the preflexion stage (41%), followed by larvae at the flexion (31%) and postflexion stages (28%, Table 1 ). Most B. robinsi larvae (99) were identified at 16 out of 52 surveyed stations from the XIXIMI-04 research cruise on oceanic water, while only 6 larvae were found at 3 out of 26 surveyed stations from the GOMECC cruise near the continental shelf ( Supplementary Table 1 ).

From the 121 samples from B. robinsi (105 larvae in different developmental stages and 16 adult tissues), 76 individuals and 11,912 loci passed the Stacks filters. After the VCFtools filters 1,140 loci and 66 individuals passed the filters. All of the scored loci were neutral (no outliers), while 97 significantly deviated from HWE and were eliminated, nine loci pairs showed significant linkage disequilibrium, and one locus of each pair was discarded. A total of 1,034 loci from 66 individuals passed ddRAD filters and were used for bioinformatic analyses: XIXIMI-04 = 55 larvae, GOMECC-3 = 3 larvae, and donated tissues = 8 adults ( Table 1 and Supplementary Table 1 ). These individuals were divided into three regions ( Figure 3 ): southern = 39, central =18, and northern= 9. The regions were based on the Miron et al. (2017) division of oceanographic provinces ( Figure 1 ). There were missing data in the final data set. In the southern region data, there were 12% missing data in 18 loci; in the central region data, there were 10% missing data in 31 loci; in the northern region data, there were 12% missing data in 282 loci, 25% missing data in 25 loci and 37% missing data in three loci. In the individuals, we registered a mean value of missing data of 1.91% (minimum = 0.19%; maximum 15.86%).

The genetic diversity descriptors obtained were very similar in the three sampled regions. Observed and expected heterozygosities showed similar values, indicating that the sampled organisms have the expected heterozygosity for a population with random mating. Also, Fis values near zero suggest that lack of endogamy in the individuals sampled within each region.

The three sampled regions showed similar allele frequencies based on non-significant Fst values. These values were consistent with the evidence of one genetic cluster observed in the PCA and DAPC graphs. These results are similar to the Fst values observed in the Western Atlantic lionfish (Pterois volitans; PLD ~25-35 days) among subpopulations located in the Northern GOM (Louisiana, Mississippi, Alabama, and Florida) (Fst = -0.040 to 0.013, p-values>0.05) suggesting the existence of one genetically homogeneous population (Pérez-Portela et al., 2018).

We did not find evidence for collective dispersal. We determined that all sampled larvae were genetically unrelated, even in stations where we sampled multiple larvae from the same developmental stage. The lack of first- and second-degree relatives and the absence of signs of inbreeding within each region where the larvae were sampled along with the existence of genetic homogeneity among all sampled regions support the lack of collective dispersal behavior. In a study of the slipper limpet (Crepidula fornicata), Riquet et al. (2017) observed a possible scenario of collective dispersal. They reported related larvae (full-sibs and half-sibs) sampled simultaneously and the presence of genetic differentiation between distinct larval groups, indicating a reduced effective population size or sweepstake reproductive success. Because we sampled unrelated larvae in different developmental stages at the same sampling stations, we inferred that these larval groups represent a mixture of larvae from different cohorts and sites, resulting from the high dispersal potential of B. robinsi consistent with a long PLD. A PLD of 4-6 months allows high rates of mixing over hundreds of kilometers, reducing the effects of sweepstakes’ reproductive success (Ottmann et al., 2016).

As consistently reflected in the results showing a lack of genetic structure (Fst, PCA, DAPC, Structure Bayesian clustering), the larvae sampled showed genetic homogeneity, indicating an open population with high degree of genetic connectivity in the GOM. Bothus robinsi larvae occurred jointly in the upper 100-meter water layer (Evseenko, 2008; Garbini et al., 2014), so their dispersal is influenced by the superficial currents ( Figure 2 ). We did not find support for the presence of high self-recruitment driven by the retention mechanisms reported for the Bay of Campeche (Sanvicente-Añorve et al., 2014; Díaz-Flores et al., 2017; Lara-Hernández et al., 2019). However, our results agree with larval dispersal simulations where the Bay of Campeche is considered a zone with high potential connectivity among nearby subpopulations (Veracruz, Tabasco, and Campeche; subpopulations located in the Campeche province and Bay of Campeche) for organisms with a PLD between 4-7 months (Sanvicente-Añorve et al., 2014; Lara-Hernández et al., 2019). These results highlight the complex dynamics of metapopulation processes within a region like the Bay of Campeche, which could be explained if export/import and self-recruitment vary along the year resulting in distinct metapopulations scenarios for different species that display different configurations of spawning times and PLDs.

The high potential connectivity results from the combination of the currents caused by the cyclonic gyre in the Bay of Campeche and the current’s back-and-forth behavior in the narrow western shelf in Veracruz and Tamaulipas ( Figure 2 ) (Sanvicente-Añorve et al., 2014; Lara-Hernández et al., 2019). The west neritic surface circulation over the narrow shelf display seasonal variations driven by the dynamics of winds, causing currents to shift towards the northwest in Spring and Summer (March to August) or southeast in Autumn and Winter (September to February) (Sanvicente-Añorve et al., 2014). This circulation keeps the larvae inside the narrow western shelf and transports larvae to the southeast (Campeche Banks) or the northwest (Flower Garden Banks). Additionally, larvae can be transported to oceanic waters towards the north through a confluence area of neritic and oceanic waters over the western edge of the Campeche Bank influenced by the cyclonic gyre of the Bay of Campeche. The gyre keeps larvae in the Bay of Campeche until their incorporation into the Campeche shelf or their possible ejection to other GOM provinces (Sanvicente-Añorve et al., 2014; Lara-Hernández et al., 2019).

For larvae with long PLD as the lobster Panulirus argus (4-7 months), the dynamic of currents in the Campeche Bay generates high values of potential connectivity. The highest connectivity values between GOM subpopulations were observed within the Campeche Bay: Campeche-Central Veracruz and Campeche-Southern Veracruz, followed by Central Veracruz-Southern Veracruz, Northern Veracruz-Central Veracruz, and Northern Veracruz-Campeche (Lara-Hernández et al., 2019). The dynamics of the Bay of Campeche, in combination with the high potential of dispersal and the multiple spawning events of B. robinsi, enhance the possibility of mixing cohorts in the dispersal, homogenizing allele frequencies within groups of larvae. Furthermore, its long PLD makes it possible for groups of larvae trapped inside the cyclonic gyre to have more possibilities to settle in the continental platform inside or outside of the Bay of Campeche (Lara-Hernández et al., 2019). The P. argus larval dispersal simulation also showed that the Bay of Campeche registered the highest values of self-recruitment (Northern Veracruz 0.711%; Central Veracruz 1.736%; Southern Veracruz 1.179%; Campeche 4.451%) (Lara-Hernández et al., 2019).

The results from studies that conducted larval dispersal simulations for P. argus (PLD=4-7 months) and the Red snapper (Lutjanus campechanus; PLD=31 days) suggest that bidirectional connectivity could be present between Florida and the Bay of Campeche for B. robinsi (PLD= approx. 5-12 months) (Evseenko, 2008; Johnson et al., 2013). In P. argus and L. campechanus larval dispersal simulations, potential connectivity is observed between all the sites located in the Bay of Campeche (Northern Veracruz, Southern Veracruz, Central Veracruz) and all the sites in Florida (Northern, Central Florida, and Southern Florida) (Lara-Hernández et al., 2019). The B. robinsi biological characteristics (subpopulation locations and PLD) are enough to have similar dispersal pathways.

We must emphasize that the direct connectivity between the Bay of Campeche and Florida is not the only way to exchange genes. The indirect exchange of genes between all provinces and subpopulations is possible under a scenario where GOM subpopulations are connected with other neighboring sites creating connectivity over multiple generations and enhancing the gene exchange via a stepping-stone model (Maruyama, 1970; Sanvicente-Añorve et al., 2014; Díaz-Flores et al., 2017; Lara-Hernández et al., 2019; Compaire et al., 2021)The dispersal pathway where the Bay of Campeche is a larval source suggests that the larvae can settle in almost any other locations within the GOM.

Another scenario that can explain genetic homogeneity is the existence of one single larval pool for the GOM. Larval dispersal simulations, along with B. robinsi long PLD and its distribution around the GOM continental shelf are consistent with the idea that B. robinsi larvae spend enough time in oceanic waters to enhance the mixing process. Distinct larval cohorts originating in the continental shelf of the Bay of Campeche, Florida, Quintana Roo, Texas, and Flower Garden Banks can contribute to the larval pools found in the central and southern regions. This mix of larvae from different areas could explain the genetic homogeneity found among our samples. We suggest that this mix of larvae is happening in all the GOM, generating the existence of one larval pool and the genetic homogenization in the GOM.

These results indicate that B. robinsi has the biological characteristic of avoiding the dynamics of retention zones. Additionally, our study confirmed the genetic homogeneity between B. robinsi larvae from the Bay of Campeche, oceanic waters from the central region of the GOM, and adults from Florida, suggesting high genetic connectivity and the existence of a single B. robinsi panmictic population.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material .

Ethical review and approval was not required for the animal study because we obtained larvae from zooplankton nets. The animals were not used in any experiment.

JM-P, AM-V, CG-S, SPJ-R, and ARO conceived and designed research. JM-P and CG-S conducted the molecular identifications. JM-P and SPJ-R conducted morphological identifications. JM-P and AM-V built the ddRAD-Seq library. Writing-original draft JM-P and AM-V. Writing-review and editing JM-P, CG-S, AM-V, SPJ-R, and ARO. JM-P and AM-V conducted the bioinformatic analysis. JM-P, CG-S, AM-V, SPJ-R, and ARO interpreted the results. All authors contributed to the article and approved the submitted version.