Climatic change drives dynamic source–sink relationships in marine species with high dispersal potential

Abstract While there is now strong evidence that many factors can shape dispersal, the mechanisms influencing connectivity patterns are species‐specific and remain largely unknown for many species with a high dispersal potential. The rock lobsters Jasus tristani and Jasus paulensis have a long pelagic larval duration (up to 20 months) and inhabit seamounts and islands in the southern Atlantic and Indian Oceans, respectively. We used a multidisciplinary approach to assess the genetic relationships between J. tristani and J. paulensis, investigate historic and contemporary gene flow, and inform fisheries management. Using 17,256 neutral single nucleotide polymorphisms we found low but significant genetic differentiation. We show that patterns of connectivity changed over time in accordance with climatic fluctuations. Historic migration estimates showed stronger connectivity from the Indian to the Atlantic Ocean (influenced by the Agulhas Leakage). In contrast, the individual‐based model coupled with contemporary migration estimates inferred from genetic data showed stronger inter‐ocean connectivity in the opposite direction from the Atlantic to the Indian Ocean driven by the Subtropical Front. We suggest that the J. tristani and J. paulensis historical distribution might have extended further north (when water temperatures were lower) resulting in larval dispersal between the ocean basis being more influenced by the Agulhas Leakage than the Subtropical Front. As water temperatures in the region increase in accordance with anthropogenic climate change, a southern shift in the distribution range of J. tristani and J. paulensis could further reduce larval transport from the Indian to the Atlantic Ocean, adding complexity to fisheries management.


| INTRODUC TI ON
One of the main objectives of fisheries management is to ensure sustainable exploitation and avoid stock depletion. Stocks generally respond predictably to fishery pressure and are able to recover under moderate levels of exploitation (Pitcher & Hart, 1982). However, an inadequate delineation of biological units and a lack of information on the source-sink dynamics could compromise optimal exploitation, replenishment, and sustainability of stocks .
Despite this, fisheries management units are frequently based on territorial jurisdictions rather than on biological information (Reiss et al., 2009). In many cases, as with the model species used in our study, Jasus tristani and Jasus paulensis, management strategies are developed without basic connectivity and genetic data for the ex- Commercial harvesting of Jasus tristani occurs around the Tristan da Cunha Archipelago (Tristan da Cunha, Inaccessible and Nightingale Islands) and Gough Island (Roscoe, 1979;Scott, 2017). The resource is currently managed using Operational Management Procedures, and total allowable catches (TACs) are set annually for each island.
Size limits and gear restrictions are also applied (Jeffs et al., 2013).
While annual catches increased gradually from 331 t in 2001 to 442 t in 2010 (Jeffs et al., 2013), CPUE values have declined since 2005 (Glass, 2014). The intensive harvesting of the J. tristani population on Vema seamount in the 1960s led to the resource being considered "fished out" within two years (Heydorn, 1969) with little sign of recovery since then (von der Heyden et al., 2007). The J. paulensis fishery at St Paul and Amsterdam Islands (Duhamel, 1980;Pruvost et al., 2015;Vranckx, 1973) is operated from Réunion Island (landing place) using one fishing vessel with associated motorized small boats. The fishery (operational since 1949/50) was overexploited in the 1970s (Pruvost et al., 2015), which led to the introduction of a number of management measures including the setting of a total allowable commercial catch (TACC), a minimum legal size limit, a closed season and the protection of egg bearing females (Jeffs et al., 2013). The fishery has now stabilized at a TACC of about 400 tonnes with an increase in CPUE noted between 2001 and 2010 (G. Duhamel, personal communication). J. tristani and J. paulensis are currently treated as separate species for the purposes of fishery management with some acknowledgment of spatial differences in age and size structure. For example, the minimum landing size varies between the Tristan da Cunha islands and is set at 70 mm carapace length at Tristan and Nightingale, which gives females a chance to spawn before being captured (Scott, 2017). However, spatial genetic structure and source-sink dynamics are not explicitly taken into account in the management. Failure to incorporate genetic information into management plans can result in local reduction of populations and in reduced overall productivity or decreasing CPUE values (Ovenden et al., 2015).
Advances in both genetic technologies and numerical modeling, and their integration into multidisciplinary approaches have improved our understanding of the complex factors that influence dispersal of marine species . Technological improvements and decreasing costs of sequencing (e.g., RAD-Seq), as well as powerful bioinformatics tools and statistical software accessibility, have resulted in an increasing use of high-resolution genomic markers (Davey et al., 2011). Markers, including single nucleotide polymorphisms (SNPs), have been applied in a range of marine animals, including nonmodel organisms, such as lobsters, and have proved to be useful for species delineation (e.g., Georges et al., 2018;Iguchi et al., 2019) and for estimating population connectivity (Silva, Villacorta-Rath, et al., 2019). Similarly, individual-based models (IBMs) have been used increasingly to infer larval dispersal and connectivity Gallego et al., 2007).
Such models are generally forced by simulated fields from oceanographic models. Increases in computational capabilities have enabled the development of more sophisticated, higher-spatial resolution oceanographic models, which in turn have allowed more detailed simulations of larval transport . Connectivity matrices from numerical models, with probabilities of larval dispersal between populations, have been combined with genetic distance and differentiation matrices in a seascape genetics approach (Selkoe et al., 2016). This multidisciplinary approach has enhanced F I G U R E 1 Sample locations and the approximate distribution range of J. tristani (blue) and J. paulensis (orange) (adapted from Booth, 2006). Triangles represent locations with both genetic data and individual-based modeling data (the latter corresponding to release sites in the individual-based model), circles represent locations with modeling data only, and the square represents a location with genetic data only our understanding of the complex drivers shaping population genetic structure (e.g., Benestan et al., 2016;Lal et al., 2017;Sandoval-Castillo et al., 2018;Singh et al., 2018;Young et al., 2015). Combining genetics with oceanographic modeling is particularly valuable for the conservation and management of marine populations as it can identify ecological and biophysical characteristics associated with particular biological units or stocks, and any asymmetries in sourcesink dynamics (Selkoe et al., 2016 and references therein).
Rock lobsters have a long pelagic larval duration (PLD; potentially up to 20 months, which has been recorded for the closely related species J. edwardsii; Bradford et al., 2015), and therefore, it has been traditionally assumed that there would be high levels of gene flow between populations. However, the concept of "open" populations is no longer considered valid as barriers to dispersal have been identified in a wide range of marine animals with medium to long pelagic larval durations (e.g., Mattingsdal et al., 2020;Silva, Villacorta-Rath, et al., 2019;Takeuchi et al., 2019;Thomas & Bell, 2013;Woodings et al., 2018). That said, species with a long PLD usually have low levels of genetic differentiation among populations as a result of large dispersal potential, large population sizes and high fecundity (Cowen & Sponaugle, 2009). Still, genetic data can also overestimate average dispersal distance as a few individuals may sporadically disperse long distances, thus reducing genetic differences among populations (Shanks, 2009). While the duration of an organism's larval phase and the extent of any environmental barriers to dispersal can be important predictors of dispersal potential, other life history characteristics, such as larval behavior and postsettlement survival, can also play a role in promoting or limiting dispersal (Cowen & Sponaugle, 2009;Moksnes et al., 2014;Villacorta-Rath et al., 2018).
Jasus tristani and Jasus paulensis distribution ranges are located over 7,000 km apart but mitochondrial DNA (mtDNA) data have shown little to no genetic differentiation between these two species (Groeneveld et al., 2012;Ovenden et al., 1997). However, the authors highlight limitations and recommend increasing sample size and the amount of mitochondrial sequence data obtained from each individual as well as including nuclear sequences. In addition, George and Kensler (1970) noted that J. tristani and J. paulensis possess a significant morphological difference in the abdominal sculpture index. Although there is mixed evidence regarding the level of the divergence between J. tristani and J. paulensis, it is clear that this species pair is genetically distinct to all other Jasus species, including J. lalandii, which is distributed on the South African coast and situated only 3,000 and 4,000 km from the J. tristani and J. paulensis populations, respectively.
The recruitment dynamics and mechanisms of connectivity between J. tristani and J. paulensis remain largely unknown as is the linkage between these two geographically distant populations. As Groeneveld et al. (2012) emphasised, the percentage of larvae dispersing from the South Atlantic to the Indian Ocean has not yet been quantified. Processes such as ocean currents (Chiswell & Booth, 1999) and larval behavior in response to environmental orientation cues (Hinojosa et al., 2016) can contribute to complex source-sink dynamics. Given the difficulty in accessing and studying these species in their natural habitat due to the isolation of the islands and challenging sea conditions, genetic and modeling approaches are valuable tools to help understand population dynamics and connectivity and inform fisheries management (Hinrichsen et al., 2011;Waples et al., 2008).
Therefore, in this study we combined genetic data obtained from high-resolution genomic markers (SNPs) with individual-based modeling that incorporates flow fields from a large-scale oceanographic model. We assessed the connectivity among J. tristani and J. paulensis populations and investigated the potential for contemporary and historical gene flow. We hypothesized that ocean currents directly influence the connectivity and genetic structure of J. tristani and J. paulensis populations. Because stocks are currently managed without any information on source-sink dynamics, results from this study have application to ongoing fisheries management decision making. Libraries were demultiplexed, and reads were filtered for overall quality (-c, -q and -r options) using process_radtags in STACKS v.2.0b9 (Catchen et al., 2013). The Stacks pipeline denovo_map.pl was executed to run each of the Stacks modules individually (ustacks, cstacks, sstacks, and populations).

| Genomic analyses
To optimize the de novo assembly, we tested a range of parameters, including m (minimum stack depth) from 3 to 7 and M (distance allowed between stacks) and n (distance allowed between catalog loci) from 1 to 9, as recommended by Rochette and Catchen (2017) and Paris et al. (2017). The parameter test showed that M = 3 and n = 3 (the M and n parameters were kept equal) provided a balance between obtaining true polymorphism and introducing sequencing error (i.e., the number of widely shared loci plateaued at about M = 3 and n = 3) and that M = 3 and n = 3 was sufficient to stabilize the proportions of loci with 1-5 SNPs. Therefore, we retained M = 3 and n = 3 for the main analysis. The high coverage with the value of m = 3 (67×) and consistent results with m = 3 to m = 7 imply a true biological signal. As m = 3 also performs well for a broad range of data sets (Paris et al., 2017;Rochette & Catchen, 2017), we retained m = 3 for the main analysis. Reads were aligned de novo with each other and a catalogue of putative RAD tags was created (cstacks module).
Putative loci were searched against the catalog (sstacks module) and further filtering was then conducted in the populations module.
Retained reads were present in at least 50% of samples within each species, detected in at least one species, had a rare allele frequency of at least 20%, and had no more than 2 alleles detected.
Potential homologs were excluded by removing markers with heterozygosity >0.50 within samples following Hohenlohe et al. (2011).
Individuals with more than 50% missing data were removed, resulting in a total of 98 samples (Table 1).
Allelic richness was estimated using hierfstat (version 0.04-22) in R (Goudet, 2005). Pairwise F ST values and respective p-values were estimated using the software GenoDive version 3.0 (Meirmans & Van Tienderen, 2004). The R package adegenet version 2.1.1 was used to estimate observed and expected heterozygosity and for estimating membership probabilities using the compoplot() function in adegenet (Jombart, 2008). Genetic structure was further investigated with principal component analyses (PCA) using the R package stats version 3.5.1 (R Core Team, 2018). Outlier analyses were used to identify signatures of selection and two complementary approaches were conducted. In BayeScan (version 2.1), prior odds were set to 100 to minimize the chance of false positives with 5,000 pilot runs, followed by 100,000 iterations (5,000 samples, a thinning interval of 10, and a burn-in of 50,000) (Foll & Gaggiotti, 2008). PCAdapt (ver-

| Migration estimates
To infer contemporary migration (first generation migrants), GeneClass2 (Piry et al., 2004) was used. In order for the GeneClass2 analysis to be computationally tractable, three subsets of 100 randomly chosen SNPs were tested (SNP subsets 1, 2 and 3). Migrants were detected by calculating the ratio of the likelihood computed from the population where the individual was sampled over the highest likelihood value among all population samples, including the population where the individual was sampled .
The assignment threshold was set to 0.01, using a Bayesian method (Rannala & Mountain, 1997). Monte-Carlo resampling method takes into account the sampling variance inherent to the limited size of the reference datasets to better control type I error rates . Monte-Carlo resampling was based on the simulation of 10,000 individuals with a Type 1 error of 0.01 over all scored loci and individuals were assigned to the location with the highest likelihood.
Historical migration (from the current generation back to the most recent common ancestor) was estimated using a coalescent method implemented in MIGRATE 3.7 (Beerli & Felsenstein, 2001).
In order for the MIGRATE analysis to be computationally tractable, three subsets of 100 randomly chosen SNPs were tested (SNP subsets 1, 2 and 3). Five independent runs were built for each subset of 100 SNPs and checked for convergence and the best run was selected. Runs were started under the SNP model using values of Θ (theta) and M (migration rate) calculated from F ST values. The first 10,000 runs were discarded as burn-in for each chain and runs consisted of sampling increments of 100 iterations with 10 replicates.
Model convergence of MIGRATE analyses was confirmed using TRACER 1.7 (Rambaut et al., 2018). The number of migrants per generation was estimated using: where Θ (theta) is a constant value per population generated from F ST calculations and M are mutation-scaled immigration rates.

| Individual-based modeling
An individual-based modeling (IBM) approach was used to investigate the transport and population connectivity of Jasus tristani and J. paulensis. The ocean flows used to drive the IBM were extracted from a global oceanographic model, specifically a 1/12° application of the Nucleus for European Modelling of the Ocean (NEMO) modeling framework, provided by the National Oceanography Centre, Southampton, UK. NEMO is a versatile modeling system that has been used for a range of modeling studies from regional to global scales, and the ability of the 1/12° configuration to represent key oceanographic features has been demonstrated (e.g., Deshayes et al., 2013;Marzocchi et al., 2015;Sérazin et al., 2015). Full details of the model may be found in Marzocchi et al. (2015). In summary,  Fox et al., 2006;Young et al., 2012). This model has been applied previously to studies of zooplankton transport (Young et al., 2014), and for the simulation of fish egg and larval transport (Fox et al., 2006;Young et al., 2015Young et al., , 2018. Although HAL includes modules allowing the simulation of temperature-dependent growth and mortality, diel vertical migration and directed swimming, there is insufficient information on the early life stages of J. tristani for inclusion of such complexity, and the model implementation was therefore comparatively simplistic. Specific model parameterizations for the simulation of Jasus sp. larvae considered spawning location and timing, depth distribution of larvae, and the duration of the pelagic larval stage. In summary, particles were advected at each model time step (30 min) according to the imposed three-dimensional velocity field, using a second-order Runge-Kutta method, with additional horizontal and vertical diffusions included using a random-walk approach (Dyke, 2001). Particles were released at six locations (JTR_Go, JTR_Tr, Ve, S1, JPA_Am, JPA_SP) with known populations of Jasus tristani and J. paulensis (Groeneveld et al., 2012; Due to limited observations of J. tristani and J. paulensis, the choice of parameterization was based on observations and prior modeling studies of J. edwardsii and J. lalandii, which have been studied in considerably more detail (e.g., Bradford et al., 2015;Bruce et al., 2007;Pollock, 1986;Stanley et al., 2015). However, the assumed spawning and planktonic periods agree with the limited observations of J. tristani phyllosoma and pueuruli (Vladimir Laptikhovsky, CEFAS, personal communication). Diel vertical migration (DVM) was not included in the representation of larval transport; there are no data on DVM in J. tristani or J. paulensis phyllosoma, and observations of DVM in J. edwardsii suggest that well-defined DVM is only observed in late-stage larvae (Bradford et al., 2005). In total, >5.7 million simulated particles were released over the 10 year study period.
For each simulation year, the percentage of larvae from each release site successfully recruiting to one of the defined recruitment zones was calculated, and the results were combined into connectivity matrices, M t , for each year describing the proportion of individuals arriving in a destination population (rows) from a given source population (columns). Summary matrices were then calculated describing the mean and standard deviation in connectivity, and the number of years with nonzero connectivity. Particles transported among J. tristani and J. paulensis sites (and vice versa) were identified and their trajectories analyzed to identify dominant inter-species connectivity pathways.

| Demographic modeling
Several hypotheses of divergence modes were tested, aiming to identify the best divergence mode of J. tristani and J. paulensis and the direction of contemporary and historical migration. Six models were built considering possible scenarios: (SI) Strict Isolation, where the environment (e.g., sea level change and ocean currents) Secondary Contact with two recent gene flow events after past isolation. All models were implemented allowing for asymmetric migration rates (m12, m21). In the models with two events of secondary contact (PAM and PSC), four parameters for migration rates were included (mA = m12 and mB = m21 in the first migration event and mC = m12 and mD = m21 in the second migration event) to determine changes in the direction of gene flow over time.
Demographic inference was performed using a folded joint site frequency spectrum (JSFS) in the diffusion approximation method implemented in the software ∂a∂i (Gutenkunst et al., 2009). The function vcf2dadi in the R package radiator (Gosselin, 2017) was used to create ∂a∂i SNP input files. Models were fitted using 20 replicate runs per model and the best model was retained. The Akaike information criterion (AIC) was used to perform comparisons among models (Sakamoto et al., 1986).
To compare among nested models of increasing complexity and address over-parameterization issues, we used the comparative framework of Tine et al. (2014) by penalizing models which contain more parameters. For each species pair, a score was estimated for each model using: where ∆ max corresponds to the difference in AIC between the worst and the best performing model (∆ max = AIC max − AIC min ) and ∆AIC i = AIC i − AIC min . Therefore, the worst model has a score of 0 and the best model has a score of 1. To evaluate the relative probabilities of the different models within each species pair, Akaike weights (W AIC ) were also calculated following: where R corresponds to the total number of models considered (R = 6).

| Migration estimates
Contemporary migration estimates showed stronger connectivity from west to east, with 71 to 86% of first generation of migrants being transported from the Atlantic (J. tristani) to the Indian Ocean (J. paulensis). The total number of first generation migrants from J. tristani to J. paulensis varied from 12 to 17 while migrants from J. paulensis to J. tristani varied from 2 to 7 across the different SNP subsets (note that 61% of migrants were commonly identified across the three SNP subsets; Table 3). Historical migration estimates showed stronger connectivity from east to west, with 82%-86% of migrants

| Individual-based modeling
The mean connectivity matrix ( Figure 3a) suggests asymmetric dispersal of Jasus larvae both between and within the ocean basins.
Predicted transport among islands in the Atlantic is strongly south- Ocean continue into the Atlantic Ocean through Agulhas leakage and shedding of Agulhas rings (Lutjeharms, 2007). Turbulent mixing in the Cape Basin region (Boebel et al., 2003) then entrains particles in the south Atlantic gyre (Figure 4b). Thus, there is a potential weak east to west stepping stone connectivity from islands in the Indian Ocean to Vema Seamount. However, the model results show stronger inter-ocean connectivity from west to east, thus suggesting that contemporary connectivity is predominantly through transport of Atlantic J. tristani larvae to Indian Ocean populations of J. paulensis.

| Demographic modeling
The

| D ISCUSS I ON
In this study, we found low but significant genetic differentiation between the rock lobsters J. tristani and J. paulensis using high-resolution genomic markers. The individual-based model and contemporary migration estimates inferred from genetic data showed stronger inter-ocean connectivity from the Atlantic to the Indian Ocean (west to east), while historical migration estimates showed stronger connectivity from east to west ( Figure 6). Results from this study have important implications for fisheries management, in particular by providing novel evidence of source-sink dynamics across ocean basins and by giving insight into the possible effects of future environmental change on population connectivity.

| Species delineation
In our study, J. tristani and J. paulensis individuals were clustered into separate groups, with low genetic differentiation and some evi-   Paul Islands have been grouped in the same zoogeographic province, called the West Wind Drift Islands Province as there is evidence of the same fish species that are endemic to both regions (Collette & Parin, 1991). Therefore, the similarities in habitat type in this zoogeographic province combined with high potential for dispersal may have facilitated gene flow between J. tristani and J. paulensis.
However, using a high number of genomic markers (17,256 SNPs), our study suggests low but significant differentiation among species.

| Historical and contemporary connectivity
In this study, historical migration estimated from the genetic data show stronger connectivity from east to west (from the Indian to the Atlantic Ocean), which is reflected in the higher genetic diversity in the Indian Ocean, and in the admixture analyses (showing introgression of J. paulensis into J. tristani). Groeneveld et al. (2012) also detected higher levels of genetic diversity in the Indian Ocean populations suggesting that these (J. paulensis) are the putative ancestral populations. However, levels of genetic diversity can also reflect higher effective population sizes as there are probably more populations of J. paulensis distributed throughout the southern Indian Ocean, whereas J. tristani may be less widely distributed in the southern Atlantic (Booth, 2006).
Compared with other lobster species (e.g., Panulirus spp. and In addition, decreasing oceanic temperatures and lower sea levels would have created newly available shallow benthic habitat for lobsters, allowing the colonization of northern areas (e.g., along the Walvis Ridge; Schaaf, 1996). Due to the high degree of uncertainty TA B L E 4 Results of model fitting for the different demographic scenarios: Strict Isolation (SI), Isolation-with-Migration (IM), Ancient Migration (AM), Secondary Contact (SC), Ancient Migration with two periods of ancient gene flow (PAM) and Secondary Contact with two periods of contact (PSC)   (Kohfeld et al., 2013). The most geographically comprehensive study (Gersonde et al., 2003) suggests that the STF may have shifted equatorward by ~3° in the Atlantic; however, such a shift relates to the movement of specific isotherms and not temperature gradients, and it is known that identifying fronts, and their associated currents, by means of a single isotherm can be misleading (Chapman et al., 2020).
Multiple proxies have been used to infer a reduction in Agulhas leakage at the LGM (e.g., Bard & Rickaby, 2009;Peeters et al., 2004).
However, Franzese et al. (2006)   landings and force fishery vessels to travel longer distances to exploit new stocks. In addition, the shift of lobsters to more southern areas could lead to marked changes in the benthic ecology and community structure in these areas (Blamey & Branch, 2012;Cockcroft et al., 2008).

| Current implications for fisheries management and considerations for the future
Ocean currents can strongly influence the genetic structure of species with long pelagic larval duration. Shifts in circulation patterns due to environmental change may affect connectivity, adding complexity to fisheries management (e.g., Potts et al., 2014). Ocean circulation is projected to change worldwide by the end of the century with changes in the intensity and position of western boundary currents including the Agulhas Current, the retroflection of which is projected to intensify and shift to the southeast (van Gennip et al., 2017;Popova et al., 2016). There is, however, considerably uncertainty regarding projected changes to circulation patterns in the Agulhas region (e.g., Tim et al., 2019), in part due to the relatively coarse spatial resolution of the presently available climate simulations. Nonetheless, an intensification of the Agulhas retroflection, in combination with a southwards shift in distribution of Jasus species, could result in reduced east to west connectivity of J. paulensis and J. tristani.
Our study provides new insights into the recruitment dynamics and connectivity of marine species with high dispersal potential and highlights the importance of incorporating multidisciplinary approaches into management decisions. Population connectivity and source-sink dynamics are currently not explicitly incorporated in the management practices of J. paulensis and J. tristani fisheries.
Accounting for such source-sink dynamics, for example by applying management measures to maintain particular levels of biomass in each population in relation to that population's likely role in supporting recruitment, may result in higher average recruitment and a more successful fishery. Centre, Southampton, UK.

CO N FLI C T O F I NTE R E S T
All authors declare that they have no competing interests.