Genetic differentiation and signatures of local adaptation revealed by RADseq for a highly dispersive mud crab Scylla olivacea (Herbst, 1796) in the Sulu Sea

Abstract Connectivity of marine populations is shaped by complex interactions between biological and physical processes across the seascape. The influence of environmental features on the genetic structure of populations has key implications for the dynamics and persistence of populations, and an understanding of spatial scales and patterns of connectivity is crucial for management and conservation. This study employed a seascape genomics approach combining larval dispersal modeling and population genomic analysis using single nucleotide polymorphisms (SNPs) obtained from RADseq to examine environmental factors influencing patterns of genetic structure and connectivity for a highly dispersive mud crab Scylla olivacea (Herbst, 1796) in the Sulu Sea. Dispersal simulations reveal widespread but asymmetric larval dispersal influenced by persistent southward and westward surface circulation features in the Sulu Sea. Despite potential for widespread dispersal across the Sulu Sea, significant genetic differentiation was detected among eight populations based on 1,655 SNPs (FST = 0.0057, p < .001) and a subset of 1,643 putatively neutral SNP markers (FST = 0.0042, p < .001). Oceanography influences genetic structure, with redundancy analysis (RDA) indicating significant contribution of asymmetric ocean currents to neutral genetic variation (Radj2 = 0.133, p = .035). Genetic structure may also reflect demographic factors, with divergent populations characterized by low effective population sizes (N e < 50). Pronounced latitudinal genetic structure was recovered for loci putatively under selection (FST = 0.2390, p < .001), significantly correlated with sea surface temperature variabilities during peak spawning months for S. olivacea (Radj2 = 0.692–0.763; p < .050), suggesting putative signatures of selection and local adaptation to thermal clines. While oceanography and dispersal ability likely shape patterns of gene flow and genetic structure of S. olivacea across the Sulu Sea, the impacts of genetic drift and natural selection influenced by sea surface temperature also appear as likely drivers of population genetic structure. This study contributes to the growing body of literature documenting population genetic structure and local adaptation for highly dispersive marine species, and provides information useful for spatial management of the fishery resource.

employed a seascape genomics approach combining larval dispersal modeling and population genomic analysis using single nucleotide polymorphisms (SNPs) obtained from RADseq to examine environmental factors influencing patterns of genetic structure and connectivity for a highly dispersive mud crab Scylla olivacea (Herbst, 1796) in the Sulu Sea. Dispersal simulations reveal widespread but asymmetric larval dispersal influenced by persistent southward and westward surface circulation features in the Sulu Sea. Despite potential for widespread dispersal across the Sulu Sea, significant genetic differentiation was detected among eight populations based on 1,655 SNPs (F ST = 0.0057, p < .001) and a subset of 1,643 putatively neutral SNP markers (F ST = 0.0042, p < .001). Oceanography influences genetic structure, with redundancy analysis (RDA) indicating significant contribution of asymmetric ocean currents to neutral genetic variation (R 2 adj = 0.133, p = .035). Genetic structure may also reflect demographic factors, with divergent populations characterized by low effective population sizes (N e < 50). Pronounced latitudinal genetic structure was recovered for loci putatively under selection (F ST = 0.2390, p < .001), significantly correlated with sea surface temperature variabilities during peak spawning months for S. olivacea (R 2 adj = 0.692-0.763; p < .050), suggesting putative signatures of selection and local adaptation to thermal clines. While oceanography and dispersal ability likely shape patterns of gene flow and genetic structure of S. olivacea across the Sulu Sea, the impacts of genetic drift and natural selection influenced by sea surface temperature also appear as likely drivers of population genetic structure. This study contributes to the growing body of literature documenting population genetic structure and local adaptation for highly dispersive marine species, and provides information useful for spatial management of the fishery resource.

| INTRODUC TI ON
Considering the spatial patterns and scales of dispersal, population connectivity has key implications for management and conservation (Moritz, 1994;Palumbi, 2003). For marine organisms, connectivity is primarily driven by complex interactions between life history characteristics and environmental conditions, which influence the dynamics and persistence of populations (Cowen & Sponaugle, 2009).
The absence of apparent physical barriers in the ocean combined with high dispersal potential of most marine organisms shaped the paradigm of genetic homogeneity in the marine environment (Hauser & Carvalho, 2008). Advances in DNA sequencing technologies now provide adequate resolution to study population genetic processes through high-throughput genotyping of single nucleotide polymorphisms (SNPs) using restriction site-associated DNA (RAD) sequencing approaches (Andrews et al., 2016;Baird et al., 2008;Davey et al., 2011) such as ddRAD (Peterson et al., 2012), 2b-RAD (Wang et al., 2012), and ezRAD . Coupling genomic approaches with biophysical modeling, which simulates and tracks larval dispersal in the marine environment (reviewed in Swearer et al., 2019), has converged into a more integrative seascape genomics approach to examine processes that shape genetic variation, whether neutral or adaptive, in the marine environment (Riginos et al., 2016;Selkoe et al., 2016).
Seascape genomics studies have improved our understanding of the environmental conditions influencing population connectivity and the spatial distribution of genetic variation in the ocean. The ability to interrogate population diversity using thousands of SNPs, and to identify loci which may be under the influence of selection versus neutral loci (Davey et al., 2011;Gagnaire et al., 2015) enables examination of adaptive divergence in response to environmental factors. There is a growing body of literature documenting genetic structure either due to neutral variation or due to adaptation to environment, at finer spatial scales than expected from species dispersal potentials, among them: ocean currents (Benestan et al., 2016;Coscia et al., 2020;Dang et al., 2019;Gilg & Hilbish, 2003;Lal et al., 2017;Paterno et al., 2017;Riginos et al., 2019;Schunter et al., 2011;Teske et al., 2016;Truelove et al., 2017;Van Wyngaarden et al., 2018;Xuereb et al., 2018), temperature (Carreras et al., 2020;Chu et al., 2014;Coscia et al., 2020;Hoey & Pinsky, 2018;Sandoval-Castillo et al., 2018;Wang et al., 2013), and salinity (Sjöqvist et al., 2015).
Mud crabs (genus: Scylla De Haan, 1833) are commercially important species with a wide distribution in mangrove areas throughout the Indo-West Pacific and other tropical and subtropical regions (Alberts-Hubatsch et al., 2015). Three mud crab species were identified in the Philippines (S. olivacea, S. tranquebarica, and S. serrata) following the description of Keenan et al. (1998), with S. olivacea (Herbst, 1796) being the most abundant (Lebata et al., 2007). While adult mud crabs exhibit limited movement (Hyland et al., 1984), larvae are thought to be highly dispersive due to their long pelagic larval duration (PLD) averaging 20-30 days (Ali et al., 2020;Jantrarotai et al., 2002;Motoh et al., 1977;Thirunavukkarasu et al., 2014), which can extend up to 75 days under suboptimal conditions of temperature and salinity (S. serrata; Baylon, 2010). Although ocean currents play a huge role in larval dispersal and settlement success (Cowen & Sponaugle, 2009), survival and development of mud crab larvae strongly depends on water temperature and salinity (Baylon, 2011;Hamasaki, 2003;Hill, 1974;Nurdiani & Zeng, 2007). In the Sulu Sea basin, S. olivacea populations are distributed along regions influenced by temporally varying environmental gradients (Oppo et al., 2003) and complex sea surface circulation such as the southward-flowing Sulu Sea throughflow from the South China Sea, the westwardflowing Bohol Sea current exiting via the Dipolog Strait, and the southern Sulu Sea gyre (Han et al., 2009;Hurlburt et al., 2011).
Oceanographic features in the Sulu Sea have been suggested to act as barriers to gene flow among populations for other taxa with relatively lower dispersal potentials than S. olivacea such as the seahorse Hippocampus spinosissimus (Lourie et al., 2005), damselfish Dascyllus aruanus (Raynal et al., 2014), and sea cucumber Holothuria scabra (Ravago-Gotanco & Kim, 2019). There is limited information, however, on the genetic structure of Philippine populations of S. olivacea, with one study reporting weak but significant genetic differentiation of Philippine populations based on microsatellite loci (Paran & Ravago-Gotanco, 2017). Moreover, there are no studies to date that explicitly examined the influence of asymmetric ocean currents and environmental heterogeneity on population connectivity and genetic structure of a highly dispersive species in the Sulu Sea.
This study examined patterns of connectivity among populations of the orange mud crab (S. olivacea) in the Sulu Sea basin. Using a seascape genomics approach, we aimed to: (a) characterize genetic structure of S. olivacea across the Sulu Sea and examine spatial patterns of genetic connectivity using SNP markers generated by RADseq; (b) examine the influence of oceanographic circulation on genetic structure and connectivity of S. olivacea populations in the Sulu Sea; and (c) examine the SNP dataset for signatures of local adaptation, which may be correlated with other environmental factors.
First, we developed a biophysical model of larval dispersal parameterized using the life history characteristics of S. olivacea, to generate realistic predictions of larval dispersal and connectivity in the Sulu Sea. We combined larval dispersal estimates with empirical genetic observations at neutral loci to determine the influence of asymmetric ocean currents on spatial patterns of connectivity. Second, we performed analyses to recover loci putatively under selection to examine signatures of local adaptation. We assessed the potential impact of environmental factors, specifically sea surface temperature and K E Y W O R D S marine connectivity, mud crab, population genomics, RAD sequencing, seascape genetics, Sulu Sea rainfall (as a proxy for salinity) on adaptive divergence of S. olivacea.
The results of this study provide valuable insights into the spatial scales of dispersal, patterns of genetic structure, and the influence of environmental and evolutionary processes on population connectivity of S. olivacea in the Sulu Sea basin, to support the development of management and conservation strategies for the fishery resource.

| Larval dispersal simulation in the Sulu Sea basin
A larval dispersal model was developed to examine connectivity of S. olivacea in the Sulu Sea. Seven larval release sites were designated, coinciding with locations where samples were collected for genetic analysis, with the exception of Coron (Table 1). Larval dispersal simulations were performed using the connectivity modeling system (CMS; Paris et al., 2013). The CMS is a model that couples  (Han et al., 2009) and the year-round spawning season of Scylla in the Philippines (Arriola, 1940). A basin-scale habitat map was included in the model using the Philippine mangrove Landsat data of Long and Giri (2011), generating 159 larval settlement nodes along the boundaries of the basin. Particle release site coordinates were adjusted up to 12 km away from the coast to simulate the offshore spawning migration reported for S. olivacea (Koolkalya et al., 2006;Moser et al., 2005). The model was configured to release fifty thousand particles from each release site weekly over a period of 1 year, for a total of 18.5 million larval particles released from seven source nodes. Released particles were parameterized as competent to settle after 20 days of passive dispersal, with a maximum duration in the water column up to 30 days based on the pelagic larval duration of S. olivacea inferred from observations on the timing of peak spawning and recruitment of juveniles to estuaries (Ali et al., 2020). To account for larval mortality, the number of particles was set to be reduced by half after 4 days of release based on the reported mortality of Scylla from larval stages zoea I to III (Jantrarotai et al., 2002;Thirunavukkarasu et al., 2014). The resulting probability estimates of larval dispersal (measured as percent settlement) were postprocessed to generate a population-by-population connectivity matrix. This was done by assigning settlement nodes to their respective islands or nearby sampling locality (~75 km radius), resulting in a reduced population source-sink dataset.

| Sample collection, species identification, and total DNA extraction
Adult Scylla olivacea (n = 146) were collected from natural mangrove habitats from 8 sites in the Sulu Sea and 2 outlier locations between 2016 and 2017 (  (Dawson et al., 1998), and stored at room temperature until analysis. Specimens were also identified using species diagnostic molecular markers following the protocol of Ma et al. (2012). DNA was extracted using the GeneJET Genomic DNA Purification Kit (Thermo Scientific), following the manufacturer's instructions with some modifications. DNA concentration was quantified using a Qubit ® fluorometer. DNA quality was assessed by agarose gel electrophoresis and measurement of absorbance ratios (A260/280) using a NanoDrop™ spectrophotometer.

| Double-digest RAD (ddRAD) sequencing
Double-digest restriction site-associated (ddRAD) libraries were prepared according to Peterson et al. (2012). DNA extracts were run on an agarose gel to check for DNA quality. Samples with low molecular weight smears were further purified using paramagnetic beads

| Read processing and SNP filtering
Sequence libraries were initially demultiplexed using the internal barcodes, while reads with low-quality scores, uncalled bases, and sequences with intact adapters were removed using the module process_radtags in STACKs v2.2 (Rochette et al., 2019). Raw read quality scores (Phred33) and adapter contamination were examined through FastQC v0.10.1 (Andrews, 2010) and MultiQC v1.7 (Ewels et al., 2016). The STACKs pipeline module denovo_map.pl was used for the construction of stacks and generation of initial catalog of putative SNPs. Stack assembly parameters such as the minimum depth of coverage required to create a stack (−m) was set to 5 (default: 3), and the maximum distance (in nucleotides) allowed between stacks (−M) and the number of mismatches allowed between sample tags when generating the catalog (−n) were increased to 4 (default: 3) to increase the SNP calling confidence and to minimize missing data (Lal et al., 2016). Additional filtering steps were performed in the populations module with the following criteria: retain only the first SNP per locus, locus must be present in all populations (−p = 10), and exclude loci that were not present in at least 50% of the individuals for each population (−r = 0.5).
Postprocessing of the SNP panel was done to exclude SNPs that were not genotyped in at least 70% of the individuals across the entire dataset. Loci with up to 30% missing data were excluded using poppr v2.8.3 (Kamvar et al., 2014), to minimize the effect of missing data on population structure inference (Reeves et al., 2016). SNP markers with minor allele frequencies (MAFs) less than 0.05 across all sites were excluded, to eliminate loci with lower power to detect genetic variability (Ardlie et al., 2002;DeWoody & DeWoody, 2005).

| Identifying non-neutral SNPs
Loci were identified as being putatively neutral or under selection using three differentiation-based (F ST ) outlier detection methods, which use different underlying models: BayeScan v2.1 (Foll & Gaggiotti, 2008), Arlequin v3.5.2.2 (Excoffier & Lischer, 2010), and OutFLANK v0.2 (Whitlock & Lotterhos, 2015). BayeScan uses a Bayesian likelihood approach to estimate the posterior probability that a locus is under selection under the assumption that allele frequencies within populations follow the multinomial-Dirichlet distribution (Feng et al., 2015;Foll & Gaggiotti, 2008). We performed the analyses on 20 pilot runs with 100,000 iterations and a burn-in of 50,000 steps. SNPs were then identified as outliers using a false discovery rate (FDR) q-value threshold of 0.05. Arlequin uses a hierarchical island model, which compares observed locus-specific F ST to the observed global F ST value using coalescent simulations (Excoffier & Lischer, 2010). We performed a total of 20,000 simulations consisting of 10 simulated groups with 100 demes per group to detect loci under selection. SNP markers with significant F ST values at the 99% confidence interval (CI) limit were considered as putative outlier loci. OutFLANK uses a maximum-likelihood approach to generate the F ST distribution of neutral loci by trimming extreme values (Whitlock & Lotterhos, 2015). We ran OutFLANK with the recommended parameters: left and right trimming fraction of 0.05, minimum heterozygosity of 0.1, and identified outlier loci using an FDR q-value threshold of 0.05. To reduce the number of false positives and adopt a conservative approach to identifying putative outlier loci, only loci identified by at least two tools were considered as outliers.
Based on the results of outlier loci analyses, we generated three SNP datasets consisting of: (a) all loci; (b) putatively neutral loci; and (c) putative outlier loci only, hereafter referred to as all, neutral, and outlier loci, respectively. For putative function annotation, consensus tags associated with candidate outlier loci were queried for sequence similarity against the NCBI nucleotide (nr/nt) collection using the BLAST algorithm BLASTN v2.6.1 (Morgulis et al., 2008;Zhang et al., 2000), with an expected threshold of 0.05.

| Population genetic structure and effective population size
The three SNP panels (all, neutral, and outlier loci) were used to examine genetic differentiation and infer connectivity among populations. We used Weir and Cockerham's F ST (1984) et al., 1992) was performed to test for population structure inferred from F ST analysis and DAPC on the three SNP datasets. AMOVA was performed using poppr v2.8.3 (Kamvar et al., 2014), with significance tested using 1,000 permutations. We examined the putatively neutral dataset for patterns of isolation by distance using the gl.ibd function in dartR, which performs a Mantel test (number of permutations = 9,999) to assess correlation between log(geographic distance) and linearized genetic distance (F ST /(1 − F ST )) matrices. Geographic distances were measured as the shortest distance over water between all pairs of sites in the Sulu Sea using the igraph package for R (Csardi & Nepusz, 2006). Effective population size was estimated for each population based on putatively neutral loci only, using the linkage disequilibrium method (N eLD ) of NeEstimator v2.01 (Do et al., 2014).

| Genetic structure and environmental factors
To examine the influence of environmental variables on genetic structure of S. olivacea in the Sulu Sea basin, we used redundancy analysis (RDA), a direct gradient analysis technique, to test for significant relationships between response and explanatory variables (Legendre et al., 2011). The Hellinger-transformed allele frequencies (Legendre & Gallagher, 2001)

| Larval dispersal estimates
Lagrangian simulations of larval dispersal using a parameterized biophysical model demonstrate ocean current-mediated connectivity among Scylla olivacea populations in the Sulu Sea basin (Figure 1a). BAT (26.2%). Self-recruitment was relatively higher for PPC (49.7%), while TWI exhibited complete self-recruitment (100%). There was no predicted settlement of particles released from western boundary sites to eastern sites. The larval dispersal model clearly reveals asymmetric transport across the Sulu Sea basin, with larval dispersal predominantly southward, and from eastern to western boundary sites.   Table S3). The divergence of CRN, PPC, and GSC was also evident for neutral loci, although F ST estimates were slightly lower (Table S4) (Table S3), but no further structure at neutral loci (Table S4).

| Genetic differentiation and population structure
The However, the discordance between pairwise F ST versus groupings recovered by DAPC and GENELAND for the southern sites BAT, NEG, and TWI populations may be influenced by missing data (19% over the outlier dataset). Thus, the clustering for these three populations should be approached with caution.
Considering the small number of outlier loci, missing data may have a big impact on the spatial genetic structure recovered by F ST and multivariate methods. To examine this further, we performed two separate analyses to handle missing data: (a) remove genotypes (individuals) with >25% missing data; and (b) impute missing data based on population frequencies as implemented in GenoDive v3.0 (Meirmans, 2020) (see Figure S2 for details). DAPC analysis of both datasets (genotypes removed and missing genotypes imputed) recovered the same four groups as the original dataset including missing data: CRN-MSJ, ANT-PPCR-ROX, BAT-NEG, and TWI ( Figure S2).

| Gene flow estimates
Coalescent simulations using Migrate-n revealed that the model representing bidirectional gene flow between East and West populations had the strongest support (Table 3) Figure S3). Migration rates between TWI and the rest of the Sulu Sea sites ranged from 28% to 75% (mean = 43%), while migration rates involving the two divergent populations PPC and CRN and other populations were generally lower (m range = 32% to 69%, mean = 37%). No statistically significant asymmetries in gene flow patterns were detected (n = 1,000 bootstraps).

| Genetic structure and environmental factors
Directional estimates of modeled larval dispersal generated seven asymmetric eigenvector maps (AEMs) representing predicted patterns of spatial genetic connectivity in the Sulu Sea.
Using redundancy analysis (RDA), backward selection of the AEM variables identified two significant predictors (AEM6 and AEM7, with p > .05) (Table 4). Together, these two AEM eigenfunctions explained 13.3% of neutral genetic variation among sites (adjusted R 2 ; p = .035). The first RDA axis (RDA1) constituted the highest proportion of genetic variation in the response data (60.9%), which is only significant at the 10% level (p = .065),  where N e values range from 139.2 to very large (infinite) at 95% CI (Table S6).
Environmental data on SST and rainfall exhibited different levels of contribution to genetic differentiation and potential latitudinal adaptation of S. olivacea in the Sulu Sea. For SST, broadly concordant results were obtained from analyses of three explanatory SST variabilities (SST mean , SST min , and SST max ), with each analysis having four independent vectors identified consisting of months mostly during the wet season (June through November, and December; Table 4).
These explanatory variables contribute to 69.2%-76.3% of the total genetic variation among sites (R 2 adj = 0.692-0.763, p < .05). RDA1 explained the highest fraction of genetic variation comprising 87.9%-88.4%, whereas RDA2 accounted for 9.8%-10.3%. Only SST mean (p = .048) and SST min (p = .049) showed statistically significant RDA1 axes, while SST max is significant at 10% level (p = .056). Rainfall data explained a lower proportion of the variation (R 2 adj = 0.656) despite having a similar number of significant variables contributing to genetic variation. The RDA model constructed using the rainfall vectors did not reveal significant correlation with genetic data (p = .089), which suggests that SST ( Figure S4 showing SST mean only) is a stronger predictor of the observed latitudinal genetic variation than rainfall.

| Genetic structure and connectivity in the Sulu Sea
The expectation of broad larval dispersal of S. olivacea based on life history features such as offshore spawning migration and a pelagic larval duration of 20-30 days (Ali et al., 2020) is consistent with larval dispersal simulations, which model that larvae have the potential to disperse widely across the Sulu Sea, a domain spanning 800 km north to south and 600 km maximum east to west. Surface circulation features modify the directionality of dispersal, which simulations show to be greater southward and eastward across the domain.
The recovery of weak yet significant genetic differentiation among F I G U R E 5 Genetic clusters of Scylla olivacea in the Sulu Sea based on GENELAND analysis of outlier loci. Posterior probability isoclines illustrate putative genetic landscapes for the Sulu Sea domain, where sites are represented by black dots, and darker colors indicate higher probabilities of membership to each of the six clusters identified across all sampling sites. Isoclines for the outgroup populations representing two genetic clusters are not shown Note: Only those variables retained by backward selection and were significant (p < .05) are included in the final model. R 2 adj represents the adjusted coefficient of determination with p-values calculated using the analysis of variance (ANOVA; permutations = 999). The proportion of constrained RDA axes were presented in columns, and values in bold and with asterisk (*) indicate significant axes at 5% and 10% level, respectively.

TA B L E 3 Log Bayes factors (LBFs) calculated from the Bezier approximation scores for four models of migration between East and West Sulu Sea populations of Scylla olivacea
S. olivacea populations across the Sulu Sea demonstrates that dispersive life history features may not necessarily lead to widespread connectivity and genetic homogeneity. This study demonstrates the greater resolution afforded by SNP loci generated from RADsequencing approaches to detect genetic differences. Using a panel of 1,655 SNPs and a reduced set of 1,643 putatively neutral SNPs, genetic structure was detected over a relatively smaller geographic area compared with a previous study based on mitochondrial DNA sequences reporting panmixia of S. olivacea from geographically disjunct sites along the western and eastern coasts of peninsular Malaysia (Strait of Malacca and South China Sea, respectively) (Rosly et al., 2013). While weak genetic differentiation was previously reported for S. olivacea based on microsatellite loci, the geographic coverage extended beyond the Sulu Sea, with a broader coverage across the Philippine archipelago (Paran & Ravago-Gotanco, 2017).
This study adds to the growing body of literature reporting significant genetic differentiation for populations of marine organisms despite the potential for broad dispersal extending beyond the spatial scales covered by genetic sampling (Hauser & Carvalho, 2008).  (Bembo et al., 1996), shore crab Carcinus aestuarii (Schiavina et al., 2014), and peacock wrasse Symphodus tinca (Carreras et al., 2017), suggesting apparent barriers to dispersal even across distances of several hundred kilometers.
Coupling genetic analysis and oceanographic modeling approaches provide additional insights into genetic structure and connectivity of populations. For S. olivacea, geographic distance does not appear to be a factor in genetic structure across the Sulu Sea.
Instead, oceanographic circulation appears to be a more significant driver of spatial patterns of dispersal and genetic structure. In particular, the predicted genetic pattern from AEM6, which identifies PPC and TWI as separate populations, is significantly correlated with the empirical allelic frequencies based on putatively neutral loci.
This pattern is broadly consistent with the genetic analyses; that is, F ST -based approaches reveal PPC as divergent from other Sulu Sea populations (Figures 2a, 3a,b). The separation of TWI, while not supported by pairwise F ST p-values after table-wide FDR adjustment, is emergent in the DAPC plot (Figure 4a). The divergence of PPC and TWI may be due to self-recruitment. TWI is modeled to have 100% self-recruitment likely due to the southern Sulu gyre, which might promote entrainment, while self-recruitment for PPC (49.7%) is relatively higher compared with the other Sulu Sea sites (22%-24%, except for BAT at 97.8%). Self-recruitment estimates are not available for CRN, but we hypothesize high rates of self-recruitment for this site considering its location in a deep embayment, which may preclude larval dispersal offshore.
Biophysical modeling reveals asymmetrical patterns of S. olivacea larval dispersal across the Sulu Sea, with greater dispersal from eastern to western boundary populations across the Sulu Sea. Genetic data, however, do not provide unequivocal support for this pattern, as approaches to infer long-term and contemporary patterns and rates of migration both point to bidirectional gene flow as the most likely scenario. Nonetheless, greater westward dispersal may still be a plausible scenario. Coalescent analyses of gene flow reveal that while the bidirectional model has the greatest likelihood, migration rates are higher from East to West populations, and westward dispersal has greater support than a model of eastward dispersal.
Moreover, given the high levels of gene flow between eastern and western populations revealed by similar allele frequency distributions exhibiting low genetic variance (F ST = 0), the hypothetical pool of migrants is expected to exhibit the same degree of genetic similarity to either West or East populations, and the accuracy of divMigrate estimates of directionality is likely to be low. Allele frequency-based methods, particularly in moderate to high gene flow systems, fare poorly at distinguishing demographically significant connections and asymmetries (Waples, 1998). Simulated data show that under high gene flow scenarios corresponding to F ST ≥ 0.005 (which are still greater than the F ST estimates for S. olivacea across the Sulu Sea), the accuracy of directionality estimates does not exceed 50% (Sundqvist et al., 2016). Thus, taking into account the oceanographic data, combined with the challenge of estimating directionality against a background of high gene flow, asymmetrical dispersal from east to west cannot be discounted. Genetic studies for other species with similar limited adult movement, but shorter pelagic larval durations than S.
olivacea, indicate limited gene flow between eastern and western boundary populations across the Sulu Sea. Genetic structure for the seahorse Hippocampus spinosissimus (Lourie et al., 2005), damselfish Dascyllus aruanus (Raynal et al., 2014), and sea cucumber Holothuria scabra (Ravago-Gotanco & Kim, 2019) attributed limited dispersal across the Sulu Sea to a combination of oceanographic circulation features such as the Sulu Sea throughflow, the geographic distance across the Sulu Sea, and the absence of stepping-stone reef habitats across the basin, as barriers to dispersal between eastern and western boundary populations. In contrast, larval dispersal simulations for three model organisms with varied dispersal potentials, a broadcast-spawning coral Acropora millepora, sea urchin Tripneustes gratilla, and a reef fish Epinephelus sp., recovered three clusters in the Sulu Sea domain (North, Central, and Southern), but did not appear to indicate restricted dispersal between eastern and western boundary populations (Pata & Yniguez, 2019).
Population allele frequencies, while largely influenced by gene flow, may also reflect demographic changes (Whitlock & McCauley, 1999). The two most divergent populations, CRN and PPC, are characterized by low estimates of effective population size (CRN N e = 24.6-26.8, PPC N e = 10.7-11.2; Table S6), indicating the possible influence of genetic drift on allele frequency of small populations, which may lead to neutral divergence (Hare et al., 2011;Waples, 2010). Genetic divergence associated with low effective population sizes have been previously reported for other marine taxa, for example, red cusk-eel Genypterus chilensis due to high fishing pressure (Córdova-Alarcón et al., 2019), and population bottlenecks for Gadus morhua (Andreev et al., 2015) and Epinephelus marginatus (Buchholz-Sørensen & Vella, 2016). The possible causes of low effective population sizes in CRN and PPC are not known. While high exploitation rates or diminished suitable habitat area may be underlying reasons for low population sizes, additional information from fishery data and habitat surveys (e.g., mangrove cover) is needed for a conclusive determination.

| Latitudinal patterns of local adaptation of Scylla olivacea in the Sulu Sea
Environmental conditions can be agents of selection shaping the genotypic composition of local populations, with environmental heterogeneity resulting in increased adaptive potential, that is, an increased average fitness of organisms in their local environment than elsewhere (Hoban et al., 2016;Sanford & Kelly, 2011). As expected, where the latitudinal thermal cline was the steepest. This significant association between latitudinal genetic structure and environmental variation during the wet season coincides with the reported peak spawning season of Scylla species in the Philippines (Arriola, 1940;Lebata et al., 2007), suggesting a biological response to environmental clines. Fine-scale genetic structure recovered by adaptive polymorphisms likely reflects the influence of temporally variable latitudinal variations in environmental variables on S. olivacea during their early life stages, despite the potential for widespread dispersal and connectivity across the Sulu Sea.
The significant association of adaptive genetic divergence with SST variability in S. olivacea reflects the influence of temperature on life history characteristics of mud crabs. Water temperature and salinity are known key factors influencing larval development, growth, and survival of mud crabs (Baylon, 2011;Hill, 1974;Nurdiani & Zeng, 2007). Variability in temperature and salinity has been reported to influence reproductive characteristics. For instance, size at maturity in mud crabs was reported to vary with latitude, with smaller size at maturity in tropical regions hypothesized to be due to faster maturation in warmer waters (Alberts-Hubatsch et al., 2015;Quinn & Kojis, 1987;Robertson & Kruger, 1994). Similar patterns of latitudinal variation in female size at maturity and fecundity were also reported for a closely related taxa, the burrowing mud crab, Helise crasa (Grapsidae) (Jones & Simons, 1981). Scylla olivacea, in particular, is known to exhibit latitudinal variation in the seasonality of peak spawning, reported to occur from July to November at latitudes between 9°N and 11°N (Koolkalya et al., 2006;Viswanathan et al., 2019), and March to September at higher latitudes (Ali et al., 2020;Ogawa et al., 2012). Latitudinal variability in temperature and salinity, through its influence on reproduction, larval survival, and development, is thus expected to play a significant role in the dynamics, genetic structure, and persistence of populations.
Patterns of genetic differentiation associated with latitudinal gradients of temperature and salinity have been reported for several marine organisms across varying spatial scales. For instance, two major latitudinal clades were recovered in the North Atlantic snail Nucella lapillus along midcoastal Maine (between 43°N and 44°N; with water and air difference reaching up to 5-10°C), in which some of the genes involved in the genetic structure were associated with heat stress tolerance (Chu et al., 2014). A pattern of population structure was also found for a high gene flow marine fish Larimichthys polyactis in the Northwest Pacific marginal seas by using an outlier locus (e.g., heat-shock protein), which is linked to local adaptation relating to seasonal variability in temperature between two regions separated by 1-2°C thermal difference between sites (Wang et al., 2013). For a marine diatom Skeletonema marinoi, a genetic break was found between the low-salinity Baltic Sea and high-salinity North Baltic Sea populations, despite the potential for migration between metapopulations based on oceanographic connectivity (Sjöqvist et al., 2015).
This study reports a pattern of latitudinal adaptive divergence associated with SST for populations of S. olivacea across the Sulu Sea.
While we recovered a small number of putative outlier loci (12 SNP loci), these were consistently identified by two outlier loci detection methods. However, as these SNPs do not map to known functional genome regions, there is no basis to generate hypotheses regarding specific genes or gene regions potentially under adaptive selection.
Further studies utilizing a larger number of SNP loci than what we were able to generate for this study (n = 1,665 loci in total) should be able to recover a proportionally greater number of putative outlier loci. Moreover, while the redundancy analysis indicates a significant association between SST and genetic variation based on outlier loci, analysis of an expanded set of environmental variables and genetic data in a more explicit genotype-environment association analysis is recommended for further studies. These are expected to contribute to greater insight into local adaptation and underlying factors influencing genetic differentiation, which is essential to understand the dynamics of populations particularly in the context of environmental variability.

| Implications for management
Understanding spatial patterns of connectivity and the environmental drivers of local adaptation of populations represents key considerations for the design of effective, resilience-based management interventions for fishery resources. For S. olivacea, basin-scale genetic differentiation was detected at both the putatively neutral and outlier loci, reflecting the influence of evolutionary (e.g., genetic drift) and environmental processes (e.g., ocean currents, temperature, and salinity) on genotypic composition of populations in the Sulu Sea. The assessment of genetic diversity and connectivity of marine populations inferred from both neutral and outlier loci provides more holistic genetic information for fisheries management of populations (Carreras et al., 2017;Gagnaire et al., 2015;Nayfa & Zenger, 2016;Sandoval-Castillo et al., 2018;Van Wyngaarden et al., 2017). In this study, we provide genetic resources (neutral and adaptive) to support the development of policy recommendations for management and conservation of S. olivacea.
From the perspective of neutral loci, S. olivacea populations in the Sulu Sea can be considered as a well-connected metapopulation, with two divergent populations (Coron and Puerto Princesa) likely influenced by some restrictions to gene flow, but also genetic drift as a consequence of small effective population sizes. Populations with low N e are particularly vulnerable to continued loss of genetic diversity and may need to be prioritized in restoration and conservation plans such as stock enhancement programs aimed at increasing yields beyond levels supported by natural recruitment (Bell et al., 2005). Stock enhancement programs initiated for the depleted mud crab S. paramamosain fishery in Japan report promising results toward increasing catch and population sizes after more than a decade of restoration efforts (Obata et al., 2006). Thus, this study recommends the development of management and conservation plans for vulnerable populations of S. olivacea, in Coron and Puerto Princesa, which are potentially facing higher rates of local extinction due to small effective population size.

Management strategies employing translocation of individuals
should also be conducted with caution, with the view to maintain localized adaptive divergence among populations. In this context, evaluation of genetic variation using outlier markers is important, to detect signatures of local adaptation. In S. olivacea, we detected a pattern of genetic structure associated with environmental gradients such as sea surface temperature. These findings are important to consider in aquaculture practices and resource management interventions that rely on translocation of individuals across geographic locations. Successful adaptation is predicted to produce genotypephenotype-environment associations, and translocation of locally adapted individuals may result in genetic environment mismatch and have significant impacts on fitness traits particularly beyond the limits of phenotypic plasticity (e.g., Kvingedal et al., 2010;Nayfa & Zenger, 2016). Thus, genetic information from this study can be used to identify sources of broodstock, which are potentially adapted to similar local environments. For example, individuals to be used in restocking the Coron population may be sourced from Mindoro, a nearby locality that is genetically similar to Coron. Likewise, the Puerto Princesa population may be restocked using individuals from an adjacent population in Roxas, Palawan, or a population across the basin, in Antique. This process may reduce outbreeding of genetically mismatched individuals that are locally adapted to different environmental conditions, which also limits the adverse effects on fitness and survival of these populations (Edmands, 2007;Edmands & Timmerman, 2003;Gharrett et al., 1999).
Overall, the results of this study can contribute to improve existing management and conservation plans for S. olivacea in the Philippines. Scylla olivacea was among the species included in a recent fisheries ordinance establishing guidelines limiting catch, trade, and transport of crablets, juvenile, and gravid individuals across the Philippines (Fisheries Administrative Order (FAO) 264s 2020; (BFAR, 2020). While not a priority species for aquaculture because of its aggressive behavior and smaller size than S. serrata, S. olivacea is the more abundant species and represents an important fishery resource that should be maintained and protected as a source of livelihood for small-scale fishers across the Philippine archipelago. It is essential to augment genetics-based approaches with other assessments of the fishery resource, to provide further insight into spatial distributions, genetic boundaries, and local adaptation in a rapidly changing marine environment, which are critical toward the design of management and conservation strategies.  Tawi-Tawi) for sample collection. We also thank Sharon Magnuson and Chris Bird (Genomics Core Lab, Texas A&M University, Corpus Christi) for performing the RAD sequencing. We thank the reviewers whose comments and insights greatly improved the manuscript. This is MSI contribution number 483.

CO N FLI C T O F I NTE R E S T
None declared.

E TH I C A L A PPROVA L
Scylla olivacea is a commercially harvested species and was not iden-