We used 1,245 genetic samples consisting of previously published (Tatarenkov et al., 2007; Tatarenkov et al., 2012; Avise and Tatarenkov, 2015; Tatarenkov et al., 2017) and unpublished data collected between 1994 and 2014 across 56 sites. From each fish, fin clips were taken and stored in either 95% ethanol or salt saturated DMSO before extracting genomic DNA using a proteinase K method described in Tatarenkov et al. (2012). For each sample, we amplified 32 unlinked microsatellite markers (Mackiewicz et al., 2006) using multiplex reactions, described in Tatarenkov et al. (2012). After alleles were separated in a GA 3100 instrument, they were scored in Genemapper 4.0 (Applied Biosystems) and binned following Tatarenkov et al. (2010).

A 10 km buffer was placed around each sampling location using the package sf (Pebesma, 2018) to match the resolution of the coarsest Ocean General Circulation Model (OGCM; ∼9 km). Overlapping buffers were grouped into populations resulting in 20 total populations from Honduras to North Central Florida. Individuals with more than 10% missing genotypic data (i.e., more than 6 missing alleles) were removed, and populations with less than 10 individuals were removed reducing the number of populations to 17 and number of individuals to 1,120 (Figure 1; Supplementary Table S1). While the genetic samples were collected across a broad time interval, previous research found samples collected from the same population over a decade apart had low genetic differentiation (F_ST = 0.023; Tatarenkov et al., 2007). Similarly, in the Florida Keys, previous work demonstrated genetic differentiation between time points within a population was low (F_ST = 0.002) and genetic differentiation was nonsignificant (Tatarenkov et al., 2012). The distance between populations over water henceforth ‘water distance’ was calculated by masking landmasses, rasterizing the marine habitat using a value of one, and calculating the least cost path between population centroids with the R package gdistance (Etten, 2017) and rnaturalearth (Massicotte and South, 2023).

To estimate oceanic connectivity, we used Connectivity Modeling System (CMS; Paris et al., 2013) to track particles spatiotemporally in a simulated ocean current environment; here, particles represent embryos or rafting adults. At its simplest configuration, CMS is a deterministic model in that particles released from the same location at the same time will follow identical trajectories. By designating a horizontal diffusivity coefficient, a random component can be added to the motion of each particle at each timestep (Okubo, 1971). The horizontal diffusivity coefficient is scaled by the timestep and grid size to account for the sub-grid turbulent processes not resolved by the OGCM (Okubo, 1971). Ocean current data were downloaded from the Hybrid Coordinate Ocean Model (HYCOM) GLBu0.08 expt_19.1 and GOMu.0.04 expt_50.1 with a uniform 0.08° and 0.04° grid size, respectively, on a 3-h output frequency (Chassignet et al., 2007). Thus, we used 36.5 (GLBu0.08 expt_19.1) and 12.3 m²s^-1 (GOMu.0.04 expt_50.1) as horizontal diffusivity coefficients within the turbulent diffusion model of CMS (equation 3 in Paris et al., 2013) following equation four in Okubo (1971). We selected 1,000 random locations (with replacement) that were an additional 2 km offshore (12 km from the nearest sampling location) along the seaward edges of the population buffers to ensure that all particles entered the simulation without being washed ashore during the first timestep. Ten particles were released from each location daily between the years 2000 and 2010 (10 particles x 1,000 locations x 17 population groupings x 365 days x 10 years = 620.5 million particles) at random times with the location of each particle retained every half hour. Particles were tracked for 63 days, which is the average time rivulus spend in the egg before successfully hatching (33 days) plus the twice the standard deviation (15 days; Earley and Marson unpublished data). A 63-day tracking period was chosen so that the simulations include 97.5% of the observed variation in successful hatching time. Because eggs or adults in flotsam that wash ashore can be exported again with ebbing tides, particles could not become stranded when encountering land masses.

From the resulting particle trajectory files, the amount of time any particle spent within the receiving population’s 10 km buffer (Time Spent In Buffer; TSIB) was calculated and summed for each source—sink combination. Some particles were only inside the population buffer at one time point, potentially due to the time interval of sampling across the simulation or ocean currents pushed the particle rapidly through the buffer. In either case, not incorporating those particles would misrepresent connectivity; therefore, any particle with only one time point was assigned a half hour of residence which was the time interval of sampling. Because oceanic connectivity is inherently asymmetric (Pop ₁ → Pop ₂ ≠ Pop ₂ → Pop ₁ ), estimates were retained as asymmetric values for directional migration analyses and made symmetric for analyses with genetic differentiation estimates. For the symmetric dataset, TSIB was summed per population pair regardless of direction (TSIB _{Pop1 → Pop2} + TSIB _{Pop2 →Pop1} ). While the summed TSIB is referred to as symmetric, the values include natural occurring asymmetries in ocean currents because the values are summed by population pair not averaged. Hence, the symmetric dataset is mathematical symmetric but still includes asymmetric processes that maybe evolutionarily relevant. However, direct connectivity via ocean currents did not exist for every population pair or in both directions; therefore, TSIB for each population pair was divided by the total TSIB across all population pairs and used as the relative probability of successful dispersal. The relative probability of successful dispersal was used to construct an undirected and directed graph (i.e., network) in the package igraph (Csardi and Nepusz, 2006) for the symmetric and asymmetric datasets, respectively; populations were the nodes and relative probability of successful dispersal were the edges. The package igraph was used to find the most direct path with highest relative dispersal probability; however, igraph finds the shortest path by summing all potential pathways between two nodes and retaining the shortest. Therefore, all edges (i.e., relative probabilities) were natural log transformed and multiplied by −1 before using igraph to find the shortest path. The transformation produces positive edge values with high relative dispersal success having smaller numbers (i.e., shorter distances). The exponential of the shortest path multiplied by −1 was divided by the number of nodes (i.e., intermediate populations) to get the per simulation relative probability of dispersal based on time. Summing transformed relative probabilities and taking the exponential of the sum multiplied by −1 is equivalent to multiplying two relative probabilities. Hence, resulting values indicate the relative probability of the most direct path scaled by the number of intermediate nodes.

All subsequent analyses were completed within R version 4.2.2 (R Core Team, 2022). Loci were tested for Hardy-Weinberg equilibrium (HWE) across the entire dataset and per population using the package pegas (Paradis, 2010). While all loci departed from HWE in at least one population, none did so in all populations; therefore, all loci were retained for analysis. Mean observed heterozygosity, expected heterozygosity, and mean allelic richness were calculated with the packages heirfstat (Goudet and Jombart, 2022) and PopGenReport (Adamack and Gruber, 2014). To evaluate the distribution of genetic variance across the range, populations were grouped into one of five regions (West Florida: TB and CC; South Florida: EG, LK, UK and FL; East Florida: SL, IR, and NS; Bahamas: SS, EI, and LB; Central America: NC, LC, TC, TA, and UH; Figure 1), and an Analysis of MOlecular Variance (AMOVA) was implemented in the package poppr (Kamvar et al., 2014) with four hierarchal levels (region, population, sampling location, individuals) and significance was tested with 999 permutations (Excoffier et al., 1992) using ade4 (Bougeard and Dray, 2018).

Population structure was evaluated to explore the spatial patterns of genetic dissimilarity using complementary methods. We performed a Discriminant Analysis of Principal Components (DAPC; Jombart et al., 2010) using the package adegenet (Jombart and Ahmed, 2011) retaining all principal component axes and evaluating K values of 1 through 60. Results from the DAPC (Supplementary Figure SA1) identified a rapid decrease in the Bayesian Information Criterion (BIC) for K = 1 through K = 5 and appeared to asymptote near K = 30; therefore, we evaluated K = 5 through K = 30 in four complementary ways: 1) TESS3 with 20 replicates and a maximum of 100,000 iterations using the package tess3r (Caye et al., 2016); 2) sNMF with 20 replicates, 100,000 iterations, and an alpha of 100 using the package LEA (Frichot and Francois, 2015); 3) STRUCTURE with 200,000 iterations, 50,000 iteration burn-in, and default parameters for three chains (Pritchard, Stephens, and Donnelly, 2000); and 4) InStruct for 100,000 iterations and 50,000 burn-in iterations for three chains (Gao, Williamson, and Bustamante, 2007). For TESS3 and sNMF, potential K values were identified by examining the minimum cross entropy value for each K. For STRUCTURE, we implemented the Evanno method (Evanno et al., 2005) and compared the loglikelihood of the data. While neither the Evanno method (Evanno et al., 2005) nor the loglikelihood of the data have been validated for use with InStruct, we implemented them in InStruct for comparison with STRUCTURE, but also evaluated the Deviance Information Criterion (DIC; Spiegelhalter et al., 2002).

Pairwise population differentiation was evaluated with G_ST (Nei and Chesser, 1983), G′_ST (Hedrick, 2005) and Jost’s D (Jost, 2008) using the package mmod (Winter 2012). R_ST (Slatkin, 1995) was also calculated between populations using the package polysat (Clark et al., 2011). G_ST (Nei and Chesser, 1983) was chosen for comparison with other metrics and previous studies; however, G′_ST (Hedrick, 2005) and Jost’s D (Jost, 2008) were specifically chosen for use with highly polymorphic microsatellite data. G_ST assumes bi-allelic loci and is based on heterozygosity (Nei and Chesser, 1983); therefore, the maximum G_ST decreases as population genetic diversity increases and systematically underestimates genetic distance (Hedrick, 1999). G′_ST attempts to account for this downward bias by standardizing G_ST (Nei and Chesser, 1983) by scaling the observed G_ST by its theoretical maximum given the mean within-population heterozygosity (Hedrick, 2005). Jost’s D estimates genetic differentiation using allelic diversity rather than heterozygosity (Jost, 2008). R_ST was chosen for comparison as it still suffers from the biases of G_ST; however, R_ST includes the mutational distance between microsattelite markers (Slatkin, 1995). We estimated directional migration with BayesAss version 3 (Wilson and Rannala, 2003) for 30 million iterations, a 10 million iteration burn-in, and thinning interval of 10,000 for five chains using 0.25, 0.085, and 0.18 for the mixing parameters for allele frequencies (Δ_A), migration rates (Δ_M), and inbreeding coefficients (Δ_F) respectively. For all Bayesian analyses (STRUCTURE, InStruct, BayesAss), convergence was assessed with Tracer v1.7.1 (Rambaut et al., 2018) and the scale reduction factor (Brooks and Gelman, 1998).

To evaluate isolation-by-distance and isolation-by-oceanography patterns, all genetic differentiation estimates were linearized (G_ST/(1—G_ST)) and water distance log10-transformed before downstream analyses (Rousset, 1997). For normality, we took the square root of relative probability of successful dispersal via ocean currents based on TSIB, which will be referred to as OC_S and OC_A for symmetric and asymmetric estimates, respectively. Mantel and partial Mantel tests were run between all linearized genetic differentiation estimates, water distance, and OC_S after controlling for water distance in the package ecodist (Goslee and Urban, 2007). To evaluate asymmetric patterns of gene flow, Mantel and partial Mantel tests were run between directional migration estimates (BayesAss), water distance, and OC_A after controlling for water distance; here, we developed and used a custom Mantel and partial Mantel function based on code from ecodist (Goslee and Urban, 2007). The custom tests do not require symmetric matrices because they transform matrices into vectors and use those vectors to derive the null distribution used to calculate the p-value. While most packages require symmetric matrices, the test does not inherently require symmetry (Mantel, 1967). Because Mantel and partial Mantel tests can be underpowered (Legendre and Fortin, 2010; Guillot and Rousset, 2013; Legendre et al., 2015; Zeller et al., 2016), we also used a Maximum-Likelihood Population Effects model (MPLE; Clarke et al., 2002) implemented in nlme (Pinheiro et al., 2023) with cormple (Pope, 2022). MLPE models account for non-independence between pairwise comparisons by incorporating the covariance between pairwise comparisons that share a common population as a random effect to account for the non-independence of pairwise data (Clarke et al., 2002). Models were run for each combination of water distance, OC (OC_S & OC_A), and their two-way interaction as predictor variables (scaled and centered) against each linearized genetic differentiation metric (e.g., G′_ST) and the fourth root of the migration rate, transformed for normality. Models then were compared using AICc following Burnham and Anderson 2004.

Mean observed heterozygosity, expected heterozygosity, and allelic richness ranged from 0 to 0.52, 0.07 to 0.69, and 1.37 to 6.13 respectively (Supplementary Table S2). While population in Central America (e.g., TA, TC) had the highest observed heterozygosity, expected heterozygosity, and allelic richness, the Florida Key (i.e., UK and LK) has comparably high levels of genetic diversity that is considerably higher that population in the Bahamas and the rest of Florida.

The AMOVA identified significant population differentiation at each level (p = 0.001) with >10% of the variance attributed to each level (Table 1). Given that the best value of K depends on the biological question (Meirmans, 2015; Janes et al., 2017), multiple K values were tested for each clustering algorithm (Supplementary Table S3; Supplementary Figure S1). STRUCTURE, using the Evanno method, identified K = 8 as the most likely number of clusters (Supplementary Table S3; Supplementary Figure S1). While the remaining cluster algorithms (i.e., InStruct, TESS3, sNMF) all identified different values of K as the most likely, they all either had inflection points (TESS3 & sNMF) or identified (InStruct with Evanno method) K∼8 as a likely number of clusters indicating that K∼8 is useful for interpreting broad scale population structure. While K = 5 was untestable by the Evanno method, the DAPC plot does indicate broadscale differentiation between Florida and Central America (Axis 1), South Florida and New Smyrna (Axis 2), and South Florida and West Florida (Axis 3; Supplementary Figure S2).

Pairwise G_ST, G′_ST, Jost’s D, and R_ST ranged from 0.03 to 0.69, 0.132 to 0.927, 0.08 to 0.817, and 0.03 to 0.66, respectively (Figure 2; Supplementary Figures S3, S4). Regardless of genetic distance metric, the Upper Keys (UK) and Lower Keys (LK) were most genetically similar (G_ST = 0.03, G′_ST = 0.13, D = 0.08, R_ST = 0.03; Figure 2; Supplementary Figures S3, S4).With G′_ST and Jost’s D, Exuma Island (EI) and Northern Caye (NC) were most genetically dissimilar (G′_ST = 0.93, D = 0.82; Figure 2; Supplementary Figure S3), while the most genetically distinct populations measured with G_ST were Exuma Island (EI) and Indian River (IR) (G_ST = 0.69; Supplementary Figure S2). With R_ST, Exuma Island (EI) and New Smyrna (NS) were the most genetically distinct (R_ST = 0.66; Supplementary Figure S4). Excluding estimated retention rates (0.67—0.96), directional migration estimated with BayesAss (Wilson and Rannala, 2003) ranged from 0.001 to 0.175 with the lowest migration rate going from Long Caye (LC) to Exuma Island (EI; 0.001) and the highest migration rate from Everglades (EG) to Lower Keys (LK; 0.175; Supplementary Figure S5). Directional migration rate confidence intervals were calculated according to Wilson and Rannala 2003 (1.96 x standard deviation). If confidence intervals for migration estimates in both directions did not overlap zero or each other, the population pair was considered to exhibit asymmetric gene flow (Supplementary Figure S5).

Mantel tests between pairwise genetic distance metrics and water distance identified significant isolation by distance with each metric (p < 0.02; Table 2), while partial Mantel tests between genetic distance and OC_S after controlling for water distance were either marginally non-significant (G_ST & G′_ST), nonsignificant (Jost’s D) with 0.08 < p < 0.15, or significant (R_ST; Table 2). After comparing models with AICc, the most likely MLPE model for each genetic distance metric included OC_S (Table 3). However, all G_ST models and the top three R_ST models were equally likely, being within 2 AICc values (Burnham and Anderson, 2004; Table 3). For both G′_ST and Jost’s D, the full model (Genetic Distance ∼ OC_S + Distance + OC_S x Distance) was best with an AICc_w > 0.99 and separated from the next best model by >10 AICc units (Table 3). For both G′_ST and Jost’s D, all variables in the best model (Distance, OC_S, and their interaction) were significant (p < 0.03); both OC_S and the interaction between OC_S and water distance were significantly negatively associated with genetic distance, while water distance was significantly positively associated with genetic distance (Table 3).

For comparison with past landscape and seascape genetic research, we developed a custom Mantel test that does not require symmetric matrices. To test our custom function (Detailed in the Mantel—Asymmetric section of the supplied code; See Data Availability Statement below), we evaluated the function on symmetric data and compared the results with ecodist. Our function and ecodist produced the same correlation values but the custom function identified OC_S as significant in the partial Mantel tests with G_ST (p = 0.016), G′_ST (p = 0.025) and R_ST (p = 0.005) but marginally non-significant with Jost’s D (p = 0.066; Table 2), while ecodist identified OC_S as marginally non-significant (G_ST: p = 0.085; G′_ST: p = 0.092), not significant (D: p = 0.146) or significant (R_ST: p = 0.042). While our custom function and ecodist calculate significance similarly, our function first converts the matrix into a vector and permutes the vector rather than matrix thereby changing, compared to existing software, the permutation procedure on which the null distribution is generated. Our custom function identified a significant negative Mantel’s r (r = −0.23; p = 0.00002) between water distance and migration rate and a nonsignificant positive association between migration rates and OC_A after controlling for water distance (p = 0.146). The full MLPE model (Migration Rate ∼ OC_A + Distance + OC_A x Distance) was the most likely model and ∼4.6 AICc units lower than the second-best model (Migration Rate = OC_A + Distance) and ∼5.4 AICc lower than the model with water distance alone with an AICc_w > 0.86 (Table 3). In the best model, water distance (p = 0.005) and the interaction between water distance and OC_A (p = 0.009) were significantly negatively associated with migration rate, while OC_A alone was nonsignificant (p = 0.479; Table 3).

Range-wide population structure was assessed using four complementary methods (STRUCTURE, InStruct, TESS3, sNMF) with 1,120 fin clips genotyped at 32 highly polymorphic microsatellite markers collected from across rivulus’ range. With coarse resolution (K = 8), methods were partially congruent, showing overlapping clusters and assignments (Supplementary Figure S6); however, there were notable discrepancies between methods (Supplementary Figure S6). The discrepancies between clustering methods indicate the importance of comparing across methodologies before interpreting ancestry results. For example, STRUCTURE clusters the southern Bahamian populations of San Salvador and Exuma Island with Indian River located on the east coast of Central Florida (Figure 1; Supplementary Figure S6). However, all other methodologies cluster Indian River with the Florida Keys (LK and UK; Supplementary Figure S6), and G_ST identified Exuma Island and Indian River to be the most genetically dissimilar (Supplementary Figure S3). Neither of these clusters were expected because neither the Florida Keys nor Southern Bahamian populations are nearest to Indian River so, interpreting STRUCTURE results without auxiliary support could lead to inappropriate inferences, which highlights the importance of conducting complementary analyses.

At K = 8, patterns of population structure and ancestry generally recapitulate previous work (Tatarenkov et al., 2017); however, our increased sample size and focus on rivulus, as opposed to the K. marmoratus species complex that includes K. hermaphroditus, simultaneously identified large blocks of shared ancestry, novel patterns of admixture (e.g., TB → CC ← FL), and distinct population clusters (e.g., NS; Figure 1). Clustering analysis generally grouped nearby populations together but there were some notable exceptions. Populations within the Florida Keys and Florida Bay (EG, LK, and UK) cluster together with limited, albeit the most, admixture from an adjacent distinct cluster containing both CC and SL, which are populations located on the west and east coast of Florida. Therefore, Florida Keys and Florida Bay populations are genetically distinct from the nearest population groups, located to the northeast and northwest. In Central America, Utila Island (Honduras), Turneffe Atoll (Belize), and Twin Cayes (Belize) formed a cluster that shows admixture with Northern Caye (Belize). Northern Caye was split between the Utila Island, Turneffe Atoll, and Twin Cayes cluster and Long Caye (Belize). Even though Long Caye is near the other Central American populations, Long Caye shares significant ancestry only with Northern Caye, and Long Caye harbors a unique genetic cluster found within an inland freshwater pond that is disconnected from seawater (Figure 1; Turko et al., 2018; Turko et al., 2022). Colonization of less saline environments can result in rapid population divergence in marine fish (Beheregaray and Sunnucks, 2001; Kelly et al., 2006) and salinity can be a strong source of selection (Velotta et al., 2022). Hence, even at regional scales, we identify a genetically distinct cluster that could potentially be maintained by local adaptation.

While our discussion is focused on broad scale clustering results (K = 8), an AMOVA suggests significant population structure between regions (Central America, Bahamas, South Florida, West Florida, East Florida), populations within regions, sampling locations within populations, individuals within sampling locations, and within individuals (Table 1). Strong population structure is common in intertidal fish (Díez-del-Molino et al., 2016; Hurt et al., 2017; Xie et al., 2022) and species occupying highly fragmented landscapes (Browne et al., 2015; Wultsch et al., 2016). We also expect strong population structure in mixed-mating species (Mori et al., 2015), those that can outcross or self-fertilize, with more structure emerging as outcrossing rates decrease (Duminil, et al., 2007; Hodgins and Yeaman, 2019). As primarily self-fertilizing hermaphrodites, a single rivulus individual or fertilized egg could colonize a new habitat fragment, resulting in extreme genetic homogeneity. If a newly founded population is isolated for an extended period due to habitat fragmentation or restricted gene flow, that population might locally adapt and outcompete immigrants, which should reduce successful gene flow from adjacent populations even if dispersal occurs (Tobler et al., 2009; Logan et al., 2016; Barbraud and Delord, 2021). Interestingly, the northernmost rivulus populations, Tampa Bay and New Smyrna, show the least mixed ancestry across the range and the most distinct clustering even at K = 8, suggesting that the populations might have arisen through rare founder events followed by limited gene flow and local adaptation. Our range-wide assessment of rivulus population structure adds more evidence to support the role of fragmentation, and mating system in driving population structure.

Previous work did not detect significant isolation by distance; however, previous work focused on local or regional scales while the current study focused on the entire species range (Tatarenkov et al., 2012; Tatarenkov et al., 2015). Both our Mantel tests and maximum likelihood population effects models (MLPE) identify distance over water as a significant predictor of genetic differentiation (Table 2; Table 3). This demonstrates that the spatial configuration of populations and distribution of landmasses are important drivers of population structure. The MLPE models also identified oceanic connectivity (OC_S) and its interaction with water distance as important predictors of genetic differentiation (Table 3). Ocean current mediated gene flow often is attributed to physical boundaries (Truelove et al., 2017; Xuereb et al., 2018); however, in this study, the distance between populations is calculated based on marine habitat thereby including physical boundaries in the distance measure. Therefore, the significance of ocean currents in this study is not due to the additional distance caused by circumventing landmasses because the additional distance is already accounted for within our distance measure. Furthermore, our network approach enabled the estimation of indirect oceanic connectivity between populations allowing for a steppingstone pattern whereby two populations may only be connected via ocean currents through intermediate populations. By dividing our estimate of indirect connectivity by the number of intermediate nodes, we standardized our connectivity estimate by the steps required to connect two populations (Supplementary Figure S6). Importantly, OC_S includes asymmetric ocean current patterns; therefore, the results suggest that unidirectional ocean current driven connectivity may impact genetic differentiation between rivulus populations. Ocean currents influence patterns of genetic differentiation in other species (Truelove et al., 2017; Huyghe and Kochzius, 2018; Xuereb et al., 2018), but these species often have pelagic larvae that can respond to oceanic stimuli to trigger settlement or vertical migration (Kelly and Palumbi, 2010; Díaz-Ferguson et al., 2012). Adult rivulus and eggs rafting on floating debris, however, cannot guide their own long-distance dispersal success, relying completely on surface currents to wash migrants ashore. Unlike previous studies that quantified the number of particles exchanged between populations (Truelove et al., 2017; Xuereb et al., 2018), we quantified the amount of time a particle spent within the population buffer. This approach enabled us to concurrently evaluate how ocean currents affect exchange of particles between populations and the impact of local ocean currents surrounding the receiving population on particle retention. While previous work demonstrates that local ocean surface currents have strong impacts on the retention of propagules by the source population (Teske et al., 2015; Feng et al., 2016), research has not explicitly evaluated the potential for local ocean surface currents surrounding the receiving population to promote dispersal by retaining immigrants. Hence, our results suggest that studies should expand from quantifying only the number of particles exchanged between two locations to include potential impacts of local surface currents that can drive the probability of successfully settling in a new location.

While asymmetric ocean currents impact the extent to which populations are genetically differentiated (Teske et al., 2015; Xuereb et al., 2018), quantifying genetic differentiation alone fails to identify asymmetries in gene flow which can impact the future spread of potentially adaptive alleles, which would then limit opportunities for genetic rescue (Matz et al., 2018; Waterhouse et al., 2018). Ocean currents are inherently directional, and the direction of currents may change over time, for instance, seasonally or over evolutionary time scales (Seigel et al., 2008; Buffoni et al., 2015). Therefore, we explicitly tested whether directional ocean current connectivity predicts directional migration rates. Our MLPE models suggested that the full model with oceanic connectivity, water distance, and their interaction was the best model even though only water distance and its interaction with ocean currents were statistically significant (Table 3). While the coefficient of the interaction term was negative (Table 3), increases in oceanic connectivity (OC_A) were positively associated with migration rates when the distance between populations over water was low to moderate (Figure 3). Our measure of ocean current connectivity did not account for the time that particles spend out to sea before arriving at the receiving population, and rivulus eggs are likely not equally viable across the duration considered. In other species, egg hatching probability decreases as eggs age (Dahlberg, 1979) and environmental conditions during egg development can impact hatching success and the speed of development (Tubbs et al., 2005; Bobe and Labbé, 2010). Rivulus has greater hatching success at lower salinities (McCain et al., 2020), and eggs exposed to air early in development can be subsequently triggered to hatch by exposure to aerial hypoxia (Wells et al., 2015). Furthermore, our biophysical model released a uniform number of particles each day across the simulation; however, rivulus exhibit reproductive seasonality that may vary across regions (Lomax et al., 2017). Hence, the timing and number of potential dispersers might show spatiotemporal variability that influences the probability of dispersal via ocean current. Therefore, even though our integration of biophysical modeling with genetic estimates of directional migration identified ocean currents as an important driver, differential egg viability and reproductive seasonality may have increased the error surrounding our oceanic connectivity estimates and highlighted potential avenues for future work focused on directional migration in marine and intertidal species.

While distance over water and its interaction with asymmetric ocean current connectivity were significant predictors of directional migration rates, migration estimates were generally small (median <0.008) with large confidence intervals (Supplementary Figure S5). Ocean surface currents are highly stochastic because local weather conditions and marine variables (Seigel et al., 2008; Buffoni et al., 2015) can change the magnitude and direction of the current across both short and long timescales (hours to millennia; Liu and Weisberg, 2005; Robinson, 1975; Sawada and Handa, 1998). Therefore, the large confidence intervals surrounding our per generation migration rates could be the product of ocean current variability along with rivulus’ short generation time (∼100 days; Cole and Noakes, 1997). Even low rates of directional migration can have large consequences for the spatial distribution of genetic variation and the spread of potentially adaptive alleles across the landscape (Rice and Papadopoulos, 2009; Sundqvist et al., 2016). Therefore, stochastic bouts of asymmetric migration may be a property of this system, resulting in periods of isolation interspersed with large influxes of migrants.

We found low rates of directional migration overall, but high migration rates (>0.15) were detected from Everglades (EG) and Fort Lauderdale (FL) to Florida Keys (UK and LK). Importantly, such migration was not reciprocated; that is, EG and FL did not receive many migrants from Florida Keys (Supplementary Figure S5). Indeed, Florida Keys received immigrants from populations across the range but, export from Florida Keys is an order of magnitude lower if it occurs at all (Supplementary Figure S5). These data indicate that islands in the Florida Keys are sinks despite being centrally located in the range. The Florida Keys harbor high genetic diversity and allelic richness relative to the rest of the range (Supplementary Table S2), with significant population structure even within that localized area (Tatarenkov et al., 2012). Migrant retention within Florida Keys may be due to a combination of ocean currents and the spatial distribution of landmasses. Surface currents from all directions converge on the Florida Keys, which could enable rafting adults or eggs to colonize (Lee et al., 1992; Weisberg et al., 2019). The Florida Keys also have landmasses distributed throughout the seascape in staggered fashion, which could facilitate retention of dispersers. Hence, local topography and more fine scale oceanographic processes are likely influencing dispersal across rivulus’ range, enabling certain locations to retain migrants more effectively while limiting their export. In this case, the Florida Keys might preserve genetic variation from previously extirpated populations, but this variation could become trapped in the Florida Keys unless ocean currents change. While not incorporated within this study, new methodologies are being developed for coastal areas that could provide more resolution into the role of local ocean current patterns on passive dispersal success in high complex intertidal systems (Gotoh and Khayyer, 2018; Roberts et al., 2019). Unfortunately, predicting how ocean currents will change with climate comes with a significant level of uncertainty (Hirschi et al., 2020; Hu et al., 2020; Bellomo et al., 2021). Even with these methodological limitations, asymmetries in migration identified potential genetic sinks. These genetic sinks could both facilitate (via elevated genetic and perhaps phenotypic variation) local adaptation to future environmental change within the sinks and interrupt the spread of potentially adaptive alleles through the range, limiting the potential for genetic rescue.

Using 10 years of biophysical simulations, we found that ocean currents are important predictors of both genetic differentiation and asymmetric gene flow, shedding light on potential drivers of rivulus’ relatively chaotic distribution of genetic variation and demonstrating that ocean currents can modulate gene flow in an intertidal species with internal fertilization and no pelagic larvae. However, the coarse resolution of current OGCM used within our biophysical models likely missed interactions between surface currents, topography, and landmass distributions that might further influence retention and export of migrants and, as a result, metapopulation structure and the spread of genetic variation across the landscape. We suggest that the Florida Keys, although central to the range, act as a genetic sink that retains migrants but does not export individuals to other populations. This could potentially be due to fine scale oceanographic processes such as eddies or tides that are currently difficult to incorporate in biophysical models. Coincidentally, new ocean modeling approaches utilizing unstructured meshes are being developed that may offer greater resolution into hydrological dynamics in coastal and intertidal systems (Sein et al., 2017). Increased resolution and incorporation of age-dependent egg viability, differential egg survival due to environmental conditions, and phenology within the biophysical model may help to further explore the roles of local ocean current fluctuations and the ecology of rivulus in mediating genetic divergence and directional gene flow between populations. Our results illuminate a complex pattern of population structure likely driven by ocean currents that modulate directional but passive dispersal of rivulus across patchy mangrove forests.