2 Materials and methods

2.1 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 Modelling System (CMS; Paris et al., 2013). The CMS is a model that couples Lagrangian-based descriptions of ocean circulation with individual-based modeling to simulate the movement of particles each with individual behaviors parameterized from known biological traits, e.g. larval duration, mortality, settlement behavior. The model was configured using contemporary oceanographic data from the 3D Global Hybrid Coordinate Ocean Model (HYCOM; Chassignet et al., 2007) with 1/25° (~4.4 km) horizontal resolution, simulating realistic sea-surface (0 -10m) currents throughout the Sulu Sea. One year of HYCOM outputs from March 2015 to February 2016 were used to run the model, covering the reversing monsoon wind forcing in the region (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 kilometers 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 one year. This resulted in 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 Scylla spp. (Jantrarotai et al., 2002, Motoh et al., 1977, Thirunavukkarasu et al., 2014). To account for larval mortality, particles were 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 post-processed 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.
 

2.2 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 to 2017 (Table 1). Species were identified using morphological characters following the description of Keenan (1998). Tissue samples were obtained from mud crab pereopods, preserved in salt saturated DMSO-EDTA (SSDE) solution  (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 manufacturer instructions with some modifications. DNA concentration was quantified using a Qubit® fluorometer. DNA quality was examined by agarose gel electrophoresis and measurement of A260/280 using a NanoDrop™ spectrophotometer.
 

2.3 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 (SPRIselect; Beckman Coulter). Approximately 150 ng DNA from individual samples were digested with MluCI and MspI and purified using AMPureXP (Beckman Coulter). Custom barcoded adapters P1 and P2 (see Peterson et al. 2012 for sequences) were then ligated to ~50 ng of DNA. The P1 adapter includes a 5bp inline unique sequence for individual barcoding. Groups of 48 samples with unique barcodes were pooled (equal volumes of each sample), purified, and size-selected using a BluePippin system (target insert size of 400 - 500 bp). Unique external indices were added to each pool by PCR amplification. PCR products were purified, fragment sizes verified, qPCR quantified, and further pooled in equimolar quantities. Sequencing of DNA libraries were performed using the Illumina NovaSeq™ 6000 Sequencing System with S4 flow cell type. Library construction and sequencing was performed at the Genomics Core Lab, Texas A&M University, Corpus Christi.
 

2.4 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), the maximum distance (in nucleotides) allowed between stacks (-M) and the number of mismatches allowed between sample tags when generating the catalogue (-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 excluding loci which were not present in at least 50% of the individuals for each population (-r = 0.5).
Post-processing 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 and DeWoody, 2005). Hardy-Weinberg equilibrium (HWE) tests at the population level were conducted using the package pegas v0.11 (Paradis, 2010), excluding loci which exhibited significant deviation from HWE (p < 0.05) after correction for false discovery rate (FDR; Benjamini and Hochberg, 1995). To limit the influence of non-independent loci, linkage disequilibrium (LD) was tested between all SNP pairs using the package genetics v1.3.8.1.2 (Warnes et al., 2019), and SNPs at strong linkage disequilibrium (r2 > 0.8) were removed (Lee et al., 2018, Tian et al., 2009). SNPs with high observed heterozygosity (Ho > 0.6) were also dropped from the dataset (e.g. (Ackiss et al., 2018, Hohenlohe et al., 2010, Van Wyngaarden et al., 2017) using pegas v0.11 (Paradis, 2010), to eliminate loci exhibiting extremely high heterozygosity resulting from false SNP calls or assembly errors (Lee et al., 2018). All R packages were run in R version 3.5.3 (R Core Team, 2019).
 

2.5 Identifying non-neutral SNPs

Loci were identified as being putatively neutral or under selection using two differentiation-based (FST) outlier detection methods which use different underlying models: BayeScan v2.1 (Foll and Gaggiotti, 2008) and Arlequin v3.5.2.2 (Excoffier and Lischer, 2010). BayeScan identifies candidate outlier loci using a Bayesian likelihood approach to estimate the posterior probability that each locus is under selection under the assumption that allele frequencies within populations follows the multinomial-Dirichlet distribution (Foll and Gaggiotti, 2008, Feng et al., 2015). Arlequin uses a hierarchical-island model which compares observed locus-specific FST to the observed global FST value using coalescent simulations (Excoffier and Lischer, 2010). These two methods were chosen for outlier analysis due to relatively lower type I (false positive) and type II (false negative) error rates compared to other outlier tests (Lotterhos and Whitlock, 2014, Narum and Hess, 2011, Vilas et al., 2012). For BayeScan, we performed an outlier test using 100,000 iterations and a burn-in of 50,000 steps in 20 pilot runs. SNPs were then identified as outliers using a false discovery rate (FDR) q-value threshold of 0.05. For Arlequin, we used 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 FST values at the 99% confidence interval (CI) limit were considered as putative outlier loci. Results of the analyses were used to generate three SNP datasets: (1) all loci; (2) putatively neutral loci only; and (3) outlier loci detected using both methods. 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).
 

2.6 Population genetic structure and effective population size

The three SNP panels (all loci,  putatively neutral loci, and outlier loci only) were used to examine genetic differentiation and infer connectivity among populations. We used Weir and Cockerham’s FST (1984) to estimate genetic differentiation over all populations (global FST) and between populations (pairwise FST), calculated using the R packages hierfstat v0.04-22 (Goudet and Jombart, 2019) and dartR v1.1.11 (Gruber et al., 2018) respectively. Significance of FST values was tested with 10,000 bootstrap replicates, and p‐values for population pairwise comparisons were adjusted for multiple tests using the false discovery rate (FDR). Pairwise FST values were visualized using heatmaps with dendrograms generated from hierarchical clustering analysis performed using the function heatmap.2 from the gplots v.3.1.0 package for R (Warnes et al., 2020).
Spatial patterns of genetic structure were examined using a discriminant analysis of principal components (DAPC; Jombart et al., 2010) implemented in the R package adegenet v2.1.3 (Jombart, 2008). DAPC was performed using sampling location as prior information. Cross-validation was performed using the function xval.DAPC with 100 replicates to determine the number of principal components to retain to avoid issues of overfitting (Jombart et al., 2010). For the outlier loci dataset, the spatially explicit clustering program GENELAND v4.0.8 (Guillot et al., 2005) was used to estimate the number of genetic clusters and infer genetic landscapes across the Sulu Sea. We used the correlated allele frequency model with the following recommended parameters: number of possible clusters (K) were initially set from 1 to 10 with 100,000 Markov Chain Monte Carlo (MCMC) iterations, thinning of 100, burn‐in of 200, maximum rate of the Poisson process fixed to 146 (N = number of individuals), maximum number of nuclei in the Poisson‐Voronoi tessellation process fixed to 438 (N multiplied by 3) and individual samples from each location were set to the same spatial coordinates. Post-processing of MCMC outputs generated a final estimate of K, which was used as the maximum number of populations in succeeding runs performed with 10 independent replicates. Individual runs were ranked, and the run having the highest average posterior probability was used to calculate individual membership coefficients, and maps of posterior probability of membership in each K cluster.
Hierarchical analysis of molecular variance (AMOVA; Excoffier et al., 1992) was performed to test for population structure inferred from FST 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 of log(geographic distance) and linearized genetic distance (FST/(1-FST)) 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 and Nepusz, 2006). Effective population size was estimated for each population based on putatively neutral loci only, using the linkage disequilibrium method (NeLD) of NeEstimator v2.01  (Do et al., 2014).
 

2.7 Genetic structure and environmental factors

We evaluated the influence of directional ocean currents, sea surface temperature (SST), and rainfall on patterns of genetic structure of S. olivacea in the Sulu Sea basin based on putatively neutral and adaptive variation.  We used redundancy analysis (RDA), a direct gradient analysis technique, to test for significant relationships between response and explanatory variables  (Legendre et al., 2011). Hellinger-transformed allele frequencies (Legendre and Gallagher, 2001) from putatively neutral and outlier SNP datasets were used as the response variable, with three environmental factors as the explanatory variables.
The effect of ocean currents on connectivity of populations was assessed using particle dispersal estimates transformed into a set of synthetic variables known as asymmetric eigenvector maps (AEMs), from which predicted spatial patterns of genetic connectivity were generated (Blanchette et al., 2008, Riginos et al., 2019, Xuereb et al., 2018). Initially, AEM eigenfunctions were generated by creating a site-by-edge (binary) matrix of connections between all pairs of sites. Weight was then attributed to each edge, and when connectivity between a given pair of sites was greater than 0 in both directions, only the direction with the highest probability of dispersal was retained. AEM eigenfunctions were generated using adespatial v0.3-7 (Dray et al., 2020) in R v3.5.3. The contributions of individual AEMs on neutral genetic variation were calculated using redundancy analysis.
The influence of SST and rainfall variabilities on genetic structure were also evaluated through redundancy analysis using the outlier dataset, to examine signatures of local adaptation along environmental gradients. For SST, high resolution (0.25o x 0.25o) monthly mean data in the Sulu Sea domain (5-14oN, 116-124oE) from 1987-2005 were downloaded from the NOAA/OAR/ESRL PSL website at https://psl.noaa.gov/ (Reynolds et al., 2007). Similarly, fine-scale monthly mean rainfall data covering the Sulu Sea basin from 1998-2014 were obtained from the Tropical Rainfall Measuring Mission database (Huffman et al., 2007). SST and rainfall data were extracted for latitudes covering the range where samples for genetic analysis were collected.
All explanatory variables (AEMs, SST, and rainfall) were individually tested using a backward selection procedure with 999 permutations using ordistep in vegan v2.5-6 (Oksanen et al., 2019), to retain the most important explanatory variables. Selected variables were included in the model, and the adjusted coefficient of determination (R2adj) was calculated. Partial RDA was performed, and the global analysis of variance was determined using the anova function in stats v3.5.3 with 999 permutations. Significance of individual RDA axes were also assessed using anova (permutations = 999) and selected environmental vectors were fitted into the ordination using envfit (permutations = 999).
 

3 Results

3.1 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). Predicted particle settlement patterns indicate predominantly southward dispersal along the western boundary. A greater proportion of particles released from ROX settled in PPC (60.1%) compared to PPC particles dispersing northward and settling in ROX (27.3%) (Figure 1b, Table S1). predominantly. Particles from ROX and PPC settled in BAT (39.3%). Conversely, very low proportions of particles released from BAT were predicted to settle northward to ROX and PPC (0.004%). A small proportion of larvae released from BAT were transported to TWI (1.8%), and the majority of BAT particles were self-recruited (97.8%). Similar patterns of greater southward dispersal are observed along the eastern boundary populations (Figure 1b), although the predicted levels of particle connectivity were lower compared to western boundary populations. Particles released from MSJ were predicted to settle in ANT (12.70%) and NEG (9.21%), particles from ANT settled in NEG (10.4%) and TWI (2.70%), and particles from NEG settled in TWI (31.06%).
Across the Sulu Sea, dispersal simulations reveal a clear pattern of westward dispersal, with larval particles released from eastern boundary sites settling in western boundary sites (Figure 1b).  In particular, particles from MSJ and ANT settled on three eastern sites: PPC (37.7% and 45.4%), ROX (11.7% and 8.10%) and BAT (5.00% and 7.80%). Particles released from NEG settled in PPC (17.6%) and 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. Clearly, the larval dispersal model reveals asymmetric transport across the Sulu Sea basin, with larval dispersal predominantly southward, and from eastern to western boundary sites.
 

3.2 SNP filtering and identification of non-neutral loci

            A total of 661,372,129 paired-end (PE) reads from 146 individual S. olivacea libraries were processed for SNP discovery and filtering using STACKs v2.2. Out of 777,762 loci, 491,192 (63.2%) were aligned with PE contigs, with an effective per-sample mean coverage of 20.8x. Following successive filtering steps (detailed in Table 2), a final dataset of 1,655 high quality, polymorphic SNP markers were recovered. BayeScan identified 12 putative outlier loci based on posterior probabilities at the 95% Bayes factor threshold, with FST values ranging from 0.0919 to 0.5308, and positive alpha values (range = 2.02 - 4.67), suggestive of diversifying selection (Foll, 2010, BayeScan v2.0 user manual). Arlequin identified 87 putative outlier loci at the 99% confidence interval limit with FST p-values values less than 0.0099. All 12 candidate outlier loci identified using BayeScan were also detected in Arlequin, thus a separate SNP panel consisting of these 12 outlier loci was generated and designated as the outlier loci dataset. Querying the contig sequences of putative outlier loci against public domain sequences using a BLAST search resulted in significant alignments of 8 RAD tags out of 12 outlier markers to known genomic regions of other crustacean species (Portunus trituberculatus, Penaeus vannamei) and fish (Chanos chanos, Danio renio). One tag had 87.76% identity match to a microsatellite region of Japanese blue swimming crab P. trituberculatus, while the rest were found to have 80%-100% identity to known and predicted genomic DNA regions (Table S2). There were no functional gene region matches for these 12 outlier loci.
 

 3.3 Genetic differentiation and population structure

All Scylla olivacea populations exhibited significant genetic differentiation based on global estimates of FST calculated using all loci (1655 SNPs; FST = 0.0070, p = 0.001) and putatively neutral loci (1643 SNPs; FST = 0.0056, p = 0.001). Excluding outgroup sites, Sulu Sea samples (n = 116 individuals; 8 sites) still exhibited significant genetic differentiation for all loci (FST = 0.0057, p = 0.001) and putatively neutral loci (FST = 0.0042, p = 0.001).  Pairwise FST values using all loci revealed the most genetically divergent populations to be CRN (FST range = 0.0021 - 0.0219), PPC (FST range = 0.0014 - 0.0187), and one outgroup GSC (FST range = 0.0014 - 0.0219) (Table S3). The divergence of CRN, PPC and GSC was also evident for neutral loci, although FST estimates were slightly lower (Table S4). Dendrograms of pairwise FST values clearly separate CRN and PPC from the rest of the Sulu Sea populations using all loci and putatively neutral loci; all p-values for pairwise comparisons were < 0.001 except for CRN-MSJ and CRN-NEG for both marker sets (Figure 3a, b). Excluding CRN and PPC, further structure is detected among the 6 Sulu Sea sites at all loci (FST = 0.0016, p = 0.002), with significant FST between TWI-ROX at all loci (Table S3), but no further structure at neutral loci (Table S4).
The DAPC using putatively neutral loci shows separation of CRN, PPC, ANT  and TWI from the other Sulu Sea populations (Figure 4a). Retaining 55 PCs based on the cross-validation results, the first three discriminant functions explained 28.0%, 26.3%, and 16.7% of the variance, respectively. The first discriminant axis reveals the separation of CRN, PPC and ANT (Figure 4b), the second discriminant axis separates TWI (Figure 4c), and the third axis further separates PPC (not shown). While pairwise FST p-values provide strong support for the separation of CRN and PPC, and to a smaller extent TWI, there is no evident support for the separation of ANT at neutral loci. An AMOVA testing a hypothesis of three genetic groups: (1) CRN, (2) PPC, and (3) the rest of the Sulu Sea populations showed significant differentiation between groups (FCT = 0.017, p = 0.002; accounting for 1.7% of the total variance).
Mantel tests show no significant relationship between geographic distance and genetic distance for  putatively neutral loci across all eight Sulu Sea sites (Mantel r = -0.086, p = 0.667). Examining eastern and western boundary populations separately, an emergent pattern of genetic distance increasing with geographic distance is observed for the eastern boundary populations (ANT, MSJ, NEG, TWI) although the relationship is not significant (Mantel r = 0.894, p = 0.083). Using the putatively neutral loci dataset, genetically divergent populations CRN, PPC and GSC were estimated to have small effective population size (Ne: CRN = 24.6-26.8, PPC = 10.7-11.2, GSC = 9.6-10.1) compared to other localities where Ne  values range from 139.2 to very large (infinite) at 95% CI (Table S6).
Outlier loci revealed pronounced genetic differentiation across the Sulu Sea (12 SNPs; FST = 0.2390, p = 0.001). A dendrogram based on pairwise FST estimates suggests four genetic clusters in the Sulu Sea: (1) CRN-MSJ; (2) ROX-PPC-ANT; (3) BAT-TWI and (4) NEG (Figure 3c). Pairwise FST among populations within each of the 3 clusters are not significant (FDR adjusted p  > 0.05), i.e. for ROX-PPC-ANT and BAT-TWI, while between-cluster comparisons are significant (FDR adjusted p < 0.05) (Table S5). The DAPC of outlier loci (11 PCs following cross-validation) suggests four genetic clusters exhibiting limited overlap in their 95% CI ellipses, and recovered the same spatial structure as pairwise FST, except that it clustered BAT and NEG, with TWI as a divergent population (Fig 4d). The first discriminant axis (58.8% of the total variance) establishes 3 groups: ROX-PPC-ANT-TWI, BAT-NEG and MSJ-CRN. The second discriminant axis (28.25% of total variance) separates TWI and establishes it as a fourth genetic group. GENELAND recovered four genetic clusters consistent with the DAPC grouping (Figure 5). AMOVA provides further support for the concordant groupings recovered by DAPC and GENELAND, with significant differentiation among the four groups (FCT = 0.209; p = 0.001) accounting for 20.9% of the total observed variance.  No further structure is detected among samples within groups (FSC = 0.0168 , p = 0.166).  However, the discordance between pairwise FST 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 FST and multivariate methods. To examine this further, we performed two separate analyses to handle missing data: (1) remove genotypes (individuals) with > 25 % missing data; (2) impute missing data based on population frequencies as implemented in GenoDive v3.0 (Meirmans 2020) (see Figure S1 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 S1). The consistent recovery of CRN-MSJ (at 12oN), ANT-PPC-ROX (at 9oN - 11oN), as genetically distinct groups from BAT-NEG and TWI (at 5oN to 8oN) by pairwise FST, DAPC and GENELAND, provides support for a pattern of latitudinal structure of Sulu Sea populations based on outlier loci.
 

3.4 Genetic structure and environmental factors

Directional estimates of modelled larval dispersal generated seven asymmetric eigenvector maps (AEMs) representing predicted patterns of spatial genetic connectivity in the Sulu Sea. Backward selection of the AEM variables identified two significant predictors (AEM6 and AEM7, with p > 0.05) (Table 4). Together, these two AEM eigenfunctions explained 13.3% of neutral genetic variation among sites (adjusted R2; p = 0.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 = 0.065), whereas RDA2 accounted for 39.1% of the total genetic variation (p = 0.292). Although AEM6 and AEM7 vectors were both selected to construct the model, individual testing of explanatory variables revealed that only AEM6 was significant (p < 0.001). The AEM6 eigenvector modelled Puerto Princesa (PPC) and Tawi-Tawi (TWI) as separate units (Figure 2b).
Environmental data on mean SST and rainfall exhibited different levels of contribution to genetic differentiation and potential latitudinal adaptation of S. olivacea in the Sulu Sea. For SST, four independent vectors were identified consisting of months mostly during the wet season (June, and August through October; Table 4) contributing 76.3% of the total genetic variation among sites (R2adj = 0.763; p = 0.041). RDA1 explained the highest fraction of genetic variation comprising 88.0% (p = 0.048), while RDA2 accounted for 10.2% (p = 0.341). Moreover, rainfall data explained a lower proportion of the variation (R2adj = 0.656) despite having a similar number of significant variables to contribute with the genetic variation. A less significant RDA model was constructed using the rainfall vectors (p = 0.089), which suggests that mean SST (Figure S2) is a stronger predictor of the observed latitudinal genetic variation than rainfall.

4 Discussion

In this study, we employed a seascape genomics approach to examine environmental factors influencing genetic structure of Scylla olivacea populations in the Sulu Sea. Putatively neutral markers revealed weak yet significant genetic differentiation significantly correlated with genetic structure predicted from particle dispersal simulations, indicating the influence of ocean currents on gene flow. Geographic distance was not a significant predictor of genetic structure. Meanwhile, outlier loci revealed a pattern of latitudinal genetic structure suggesting local adaptation to latitudinal environmental gradients. Mean SST is a stronger predictor of adaptive divergence along environmental gradients than rainfall. The results presented here provided evidence of basin-scale genetic differentiation of S. olivacea populations in the Sulu Sea, which may be used to devise spatially explicit management and conservation interventions.
 

4.1 Genetic structure and connectivity in the Sulu Sea

Broad dispersal of S. olivacea may be inferred from life history features such as offshore spawning migration and a relatively long pelagic larval duration which may reach up to 75 days. Biophysical modeling predicts that larvae have the potential to disperse widely across the Sulu Sea, an area 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. However, the recovery of weak yet significant genetic differentiation among 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 RAD-sequencing approaches to detect genetic differences. Using a panel of 1,655 SNPs and a reduced set of 1,643 putatively neutral SNPs excluding outlier loci, genetic structure was detected over a relatively smaller geographic area compared to 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 reported for S. olivacea based on microsatellite loci, geographical coverage was broader across the Philippine archipelago and was not limited to the Sulu Sea (Paran and 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  (Hauser and Carvalho, 2008). Weak yet significant genetic differentiation, with comparable estimates of low FST values over sampling scales of hundreds of kilometers have been reported for broadly-dispersing species such as highly mobile fish (Atlantic cod, FST = 0.004; Knutsen et al., 2003), or invertebrates with extensive larval duration periods such as red rock lobsters (FST = 0.004; Iacchei et al., 2013) and spiny lobsters (FST = 0.0016; Truelove et al., 2017). In the Adriatic Sea (800 km long, 200 km wide), a semi-enclosed ocean basin comparable in area to the Sulu Sea (800 km long, 600 km wide), significant genetic differentiation was also reported for a range of organisms with similar bipartite life histories and broad dispersal potentials, such as the anchovy Engraulis encrasicolus (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. 
Genetic differentiation across the Sulu Sea does not appear to be a function of geographic distance alone, but is influenced by oceanographic features based on the significant contributions of AEM predictors on the neutral genetic variation among Sulu Sea populations (R2adj = 0.133, p = 0.035). In particular, the predicted genetic pattern from AEM6, which identifies PPC and TWI as separate populations, is significantly correlated to the empirical allelic frequencies based on putatively neutral loci. This pattern is broadly consistent with the genetic analyses, i.e. FST-based approaches reveal PPC as divergent from other Sulu Sea populations (Figure 2a, 3a, 3b). The separation of TWI, while not supported by pairwise FST p-values after table-wide FDR adjustment, is emergent in the DAPC plot (Figure 4d). 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 entraintment, while self-recruitment for PPC (49.7%) is relatively higher compared to 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 considering the deep embayment of CRN, which may preclude larval dispersal offshore. Excluding CRN, PPC and TWI, the rest of the Sulu Sea populations do not exhibit further genetic differentiation (FST = 0.0004, p = 0.094), suggesting no apparent barriers to S. olivacea larval dispersal across the Sulu Sea. 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. 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 and 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 and Yniguez, 2019).­­
Population allele frequencies, while largely influenced by gene flow, may also reflect demographic changes (Whitlock and McCauley, 1999). The two most divergent populations, CRN and PPC are characterized by low estimates of effective population size (CRN Ne = 24.6-26.8, PPC Ne = 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, e.g.  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 and 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) are needed to make a more conclusive determination.
 

4.2 Latitudinal adaptation of S. 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 i.e. an increased average fitness of organisms in their local environment than elsewhere (Hoban et al., 2016, Sanford and Kelly, 2011). Scylla olivacea populations exhibit pronounced genetic differentiation at outlier loci, demonstrating fine-scale latitudinal genetic structure across the Sulu Sea. Four genetic clusters were identified using multiple genetic approaches (FST, DAPC, GENELAND), with AMOVA indicating significant differentiation among groups accounting for 20.9% of the total variance (FCT = 0.209; p = 0.001): (1) CRN-MSJ, (2) ROX-PPC-ANT, (3) BAT-NEG, and (4) TWI. Two environmental variables, monthly mean SST and rainfall can explain the latitudinal genetic structure of S. olivacea, with SST (R2adj = 0.763, p = 0.041) as a stronger predictor of genetic variation than rainfall (R2adj = 0.656, p = 0.089). For both SST and rainfall data, most of the selected variables (months) that were included in the model covers the wet season (June - November) where the latitudinal thermal cline was steepest. Moreover, the 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 connectivity.
Temperature and salinity play a significant role in the seasonal occurrence and abundance of mudcrabs. Changes brought about by increased freshwater flow during the rainy season are thought to influence mud crab abundance through enhanced productivity of coastal areas and estuaries (Alberts-Hubatsch et al., 2015, Butcher et al., 2002, Meynecke et al., 2008). Moreover, latitudinal variations in temperature and salinity can be expected to drive variability in reproductive characteristics.  Size-at-maturity in mud crabs has been shown 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 and Kojis, 1987, Robertson and Kruger, 1994). Latitudinal variation in reproductive characteristics was also reported for a closely related taxa, the burrowing mud crab, Helise crasa (Grapsidae) where the maximum crab size, size of maturity of females, and numbers of eggs carried per female increased significantly with increased latitude (Jones and Simons, 1981).  Latitudinal variation in peak spawning of S. olivacea, which occurs from July to November at latitudes between 9oN to 11oN (Koolkalya et al., 2006, Viswanathan et al., 2019), and from March to September at higher latitudes (Ali et al., 2020, Ogawa et al., 2012)  suggests a potential for adaptive variation across broad latitudinal scales. Moreover, water temperature and salinity are known key factors influencing larval development, growth and survival of mud crabs (Baylon, 2011, Hill, 1974, Nurdiani and Zeng, 2007). Thus, latitudinal variability in temperature and salinity are expected to significantly impact reproduction, larval survival, and development, and ultimately, 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, 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 is the first report of latitudinal adaptive divergence for the mudcrab. Assessment of local adaptation and underlying factors influencing genetic differentiation is essential to understand the dynamics of populations associated with environmental factors.

4.3 Implications to Management

Understanding of spatial patterns of connectivity and 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, 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 (Sandoval-Castillo et al., 2018, Carreras et al., 2017, Gagnaire et al., 2015, Nayfa and Zenger, 2016, 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 metapopulation with two divergent populations (Coron and Puerto Princesa) likely influenced by genetic drift as a consequence of small effective population sizes. Populations with low Ne 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 towards 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 temperature and salinity. 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 genotype-phenotype-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 and 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 which 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 and 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) 264 s 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 towards the design of management and conservation strategies.

Literature cited

Ackiss, A. S., Bird, C. E., Akita, Y., Santos, M. D., Tachihara, K. & Carpenter, K. E. 2018. Genetic patterns in peripheral marine populations of the fusilier fish Caesio cuning within the Kuroshio Current. Ecol Evol, 8, 11875-11886.
Alberts-Hubatsch, H., Lee, S. Y., Meynecke, J.-O., Diele, K., Nordhaus, I. & Wolff, M. 2015. Life-history, movement, and habitat use of Scylla serrata (Decapoda, Portunidae): current knowledge and future challenges. Hydrobiologia, 763, 5-21.
Ali, M. Y., Hossain, M. B., Sana, S., Rouf, M. A., Yasmin, S. & Sarower, M. G. 2020. Identifying peak breeding season and estimating size at first maturity of mud crab (Scylla olivacea) from a coastal region of Bangladesh. Heliyon, 6, e04318.
Andreev, V., Fokin, M., Mugue, N. & Strelkov, P. 2015. Long-term persistence and evolutionary divergence of a marine fish population with a very small effective population size (Kildin cod Gadus morhua kildinensis). Marine Biology, 162, 979-992.
Andrews, K. R., Good, J. M., Miller, M. R., Luikart, G. & Hohenlohe, P. A. 2016. Harnessing the power of RADseq for ecological and evolutionary genomics. Nature Reviews Genetics, 17, 81.
Andrews, S. 2010. FastQC: a quality control tool for high throughput sequence data. Babraham Bioinformatics, Babraham Institute, Cambridge, United Kingdom.
Ardlie, K. G., Lunetta, K. L. & Seielstad, M. 2002. Testing for population subdivision and association in four case-control studies. Am J Hum Genet, 71, 304-11.
Arriola, F. J. 1940. A Preliminary Study of the Life History of Scylla Serrata (Forskal), Philippine Journal of Science.
Baird, N. A., Etter, P. D., Atwood, T. S., Currey, M. C., Shiver, A. L., Lewis, Z. A., Selker, E. U., Cresko, W. A. & Johnson, E. A. 2008. Rapid SNP discovery and genetic mapping using sequenced RAD markers. PloS one, 3, e3376.
Baylon, J. C. 2010. Effects of Salinity and Temperature on Survival and Development of Larvae and Juveniles of the Mud Crab, Scylla serrata (Crustacea: Decapoda: Portunidae). Journal of the World Aquaculture Society, 41, 858-873.
Baylon, J. C. 2011. Survival and development of larvae and juveniles of the mud crab Scylla olivacea Forskal (Crustacea: Decapoda: Portunidae) at various temperatures and salinities. The Philippine Agricultural Scientist, 94.
Bell, J. D., Munro, J. L., Nash, W. J., Rothlisberg, P. C., Loneragan, N. R., Ward, R. D. & Andrew, N. L. B. T. a. I. M. B. 2005. Restocking and stock enhancement of marine invertebrate fisheries. Academic Press.
Bembo, D. G., Carvalho, G. R., Cingolani, N., Arneri, E., Giannetti, G. & Pitcher, T. J. 1996. Allozymic and morphometric evidence for two stocks of the European anchovy Engraulis encrasicolus in Adriatic waters. Marine Biology, 126, 529-538.
Benestan, L., Quinn, B. K., Maaroufi, H., Laporte, M., Clark, F. K., Greenwood, S. J., Rochette, R. & Bernatchez, L. 2016. Seascape genomics provides evidence for thermal adaptation and current-mediated population structure in American lobster (Homarus americanus). Molecular ecology, 25, 5073-5092.
Benjamini, Y. & Hochberg, Y. 1995. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society. Series B (Methodological), 57, 289-300.
Bfar 2020. Regulation on the catching, possession, transporting, selling, trading and exporting of Mangrove Crablets, Juvenile Mangrove Crabs, and Gravid Mangrove Crabs (Scylla Spp.). In: RESOURCES, B. O. F. A. A. (ed.) FAO 264 s 2020.
Blanchette, C. A., Melissa Miner, C., Raimondi, P. T., Lohse, D., Heady, K. E. K. & Broitman, B. R. 2008. Biogeographical patterns of rocky intertidal communities along the Pacific coast of North America. Journal of Biogeography, 35, 1593-1607.
Buchholz-Sørensen, M. & Vella, A. 2016. Population Structure, Genetic Diversity, Effective Population Size, Demographic History and Regional Connectivity Patterns of the Endangered Dusky Grouper, Epinephelus marginatus (Teleostei: Serranidae), within Malta’s Fisheries Management Zone. PLOS ONE, 11, e0159864-e0159864.
Butcher, P., Boulton, A. J. & Smith, S. 2002. Mud crab (Scylla serrata: Portunidae) populations as indicators of the effectiveness of estuarine marine protected areas. World Congress on Aquatic Protected Areas Proceedings, 421-427.
Carreras, C., Ordonez, V., Zane, L., Kruschel, C., Nasto, I., Macpherson, E. & Pascual, M. 2017. Population genomics of an endemic Mediterranean fish: differentiation by fine scale dispersal and adaptation. Sci Rep, 7, 43417.
Chassignet, E. P., Hurlburt, H. E., Smedstad, O. M., Halliwell, G. R., Hogan, P. J., Wallcraft, A. J., Baraille, R. & Bleck, R. 2007. The HYCOM (HYbrid Coordinate Ocean Model) data assimilative system. Journal of Marine Systems, 65, 60-83.
Chu, N. D., Kaluziak, S. T., Trussell, G. C. & Vollmer, S. V. 2014. Phylogenomic analyses reveal latitudinal population structure and polymorphisms in heat stress genes in the North Atlantic snail Nucella lapillus. Mol Ecol, 23, 1863-73.
Córdova-Alarcón, V. R., Araneda, C., Jilberto, F., Magnolfi, P., Toledo, M. I. & Lam, N. 2019. Genetic Diversity and Population Structure of Genypterus chilensis, a Commercial Benthic Marine Species of the South Pacific. Frontiers in Marine Science, 6, 748-748.
Cowen, R. K. & Sponaugle, S. 2009. Larval dispersal and marine population connectivity. Ann Rev Mar Sci, 1, 443-66.
Csardi, G. & Nepusz, T. 2006. The igraph software package for complex network research. InterJournal, complex systems, 1695, 1-9.
Davey, J. W., Hohenlohe, P. A., Etter, P. D., Boone, J. Q., Catchen, J. M. & Blaxter, M. L. 2011. Genome-wide genetic marker discovery and genotyping using next-generation sequencing. Nature Reviews Genetics, 12, 499-510.
Dawson, M. N., Raskoff, K. A. & Jacobs, D. K. 1998. Field preservation of marine invertebrate tissue for DNA analyses. Molecular marine biology and biotechnology, 7, 145-152.
Dewoody, Y. D. & Dewoody, J. A. 2005. On the Estimation of Genome-wide Heterozygosity Using Molecular Markers. Journal of Heredity, 96, 85-88.
Do, C., Waples, R. S., Peel, D., Macbeth, G. M., Tillett, B. J. & Ovenden, J. R. 2014. NeEstimator v2: re-implementation of software for the estimation of contemporary effective population size (Ne ) from genetic data. Mol Ecol Resour, 14, 209-14.
Dray, S., Bauman, D., Blanchet, G., Borcard, D., Clappe, S., Guenard, G., Jombart, T., Larocque, G., Legendre, P. & Madi, N. 2020. adespatial: Multivariate multiscale spatial analysis. R package version 0.3-7. Available in: https://CRAN. R-project. org/package= adespatial.
Edmands, S. 2007. Between a rock and a hard place: evaluating the relative risks of inbreeding and outbreeding for conservation and management. Mol Ecol, 16, 463-75.
Edmands, S. & Timmerman, C. C. 2003. Modeling Factors Affecting the Severity of Outbreeding Depression. Conservation Biology, 17, 883-892.
Ewels, P., Magnusson, M., Lundin, S. & Kaller, M. 2016. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics, 32, 3047-8.
Excoffier, L. & Lischer, H. E. 2010. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour, 10, 564-7.
Excoffier, L., Smouse, P. E. & Quattro, J. M. 1992. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics, 131, 479-491.
Feng, X. J., Jiang, G. F. & Fan, Z. 2015. Identification of outliers in a genomic scan for selection along environmental gradients in the bamboo locust, Ceracris kiangsu. Sci Rep, 5, 13758.
Foll, M. & Gaggiotti, O. 2008. A Genome-Scan Method to Identify Selected Loci Appropriate for Both Dominant and Codominant Markers: A Bayesian Perspective. Genetics, 180, 977-993.
Gagnaire, P. A., Broquet, T., Aurelle, D., Viard, F., Souissi, A., Bonhomme, F., Arnaud-Haond, S. & Bierne, N. 2015. Using neutral, selected, and hitchhiker loci to assess connectivity of marine populations in the genomic era. Evolutionary Applications, 8, 769-786.
Gharrett, A. J., Smoker, W. W., Reisenbichler, R. R. & Taylor, S. G. 1999. Outbreeding depression in hybrids between odd- and even-broodyear pink salmon. Aquaculture, 173, 117-129.
Gilg, M. R. & Hilbish, T. J. 2003. The Geography of Marine Larval Dispersal: Coupling Genetics with Fine-Scale Physical Oceanography. Ecology, 84, 2989-2998.
Goudet, J. & Jombart, T. 2019. hierfstat: Estimation and Tests of Hierarchical F-Statistics. R package version 0.04–22. 2015.
Gruber, B., Unmack, P. J., Berry, O. F. & Georges, A. 2018. dartr: An r package to facilitate analysis of SNP data generated from reduced representation genome sequencing. Mol Ecol Resour, 18, 691-699.
Guillot, G., Mortier, F. & Estoup, A. 2005. Geneland: a computer package for landscape genetics. Molecular Ecology Notes, 5, 712-715.
Hamasaki, K. 2003. Effects of temperature on the egg incubation period, survival and developmental period of larvae of the mud crab Scylla serrata (Forskål) (Brachyura: Portunidae) reared in the laboratory. Aquaculture, 219, 561-572.
Han, W., Moore, A. M., Levin, J., Zhang, B., Arango, H. G., Curchitser, E., Di Lorenzo, E., Gordon, A. L. & Lin, J. 2009. Seasonal surface ocean circulation and dynamics in the Philippine Archipelago region during 2004–2008. Dynamics of Atmospheres and Oceans, 47, 114-137.
Hare, M. P., Nunney, L., Schwartz, M. K., Ruzzante, D. E., Burford, M., Waples, R. S., Ruegg, K. & Palstra, F. 2011. Understanding and estimating effective population size for practical application in marine species management. Conserv Biol, 25, 438-49.
Hauser, L. & Carvalho, G. R. 2008. Paradigm shifts in marine fisheries genetics: ugly hypotheses slain by beautiful facts. Fish and Fisheries, 9, 333-362.
Hill, B. J. 1974. Salinity and temperature tolerance of zoeae of the portunid crab Scylla serrata. Marine Biology, 25, 21-24.
Hoban, S., Kelley, J. L., Lotterhos, K. E., Antolin, M. F., Bradburd, G., Lowry, D. B., Poss, M. L., Reed, L. K., Storfer, A. & Whitlock, M. C. 2016. Finding the Genomic Basis of Local Adaptation: Pitfalls, Practical Solutions, and Future Directions. Am Nat, 188, 379-97.
Hohenlohe, P. A., Bassham, S., Etter, P. D., Stiffler, N., Johnson, E. A. & Cresko, W. A. 2010. Population genomics of parallel adaptation in threespine stickleback using sequenced RAD tags. PLoS Genet, 6, e1000862.
Huffman, G. J., Bolvin, D. T., Nelkin, E. J., Wolff, D. B., Adler, R. F., Gu, G., Hong, Y., Bowman, K. P. & Stocker, E. F. 2007. The TRMM multisatellite precipitation analysis (TMPA): Quasi-global, multiyear, combined-sensor precipitation estimates at fine scales. Journal of hydrometeorology, 8, 38-55.
Hurlburt, H. E., Metzger, J. E., Sprintall, J., Riedlinger, S. N., Arnone, R. A., Shinoda, T. & Xu, X. 2011. Circulation in the Philippine Archipelago simulated by 1/12° and 1/25° Global HYCOM and EAS NCOM. Oceanography, 24, 28-47.
Hyland, S., Hill, B. & Lee, C. 1984. Movement within and between different habitats by the portunid crab Scylla serrata. Marine Biology, 80, 57-61.
Iacchei, M., Ben-Horin, T., Selkoe, K. A., Bird, C. E., García-Rodríguez, F. J. & Toonen, R. J. 2013. Combined analyses of kinship and FST suggest potential drivers of chaotic genetic patchiness in high gene-flow populations. Molecular Ecology, 22, 3476-3494.
Jantrarotai, P., Taweechuer, K. & Pripanapong, S. 2002. Salinity Levels on Survival Rate and Development of Mud Crab ( Scylla olivacea ) from Zoea to Megalopa and from Megalopa to Crab Stage. Kasetsart Journal (Natural Science), 36, 278-284.
Jombart, T. 2008. adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics, 24, 1403-5.
Jombart, T., Devillard, S. & Balloux, F. 2010. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet, 11, 94.
Jones, M. B. & Simons, M. J. 1981. Habitat preferences of two estuarine burrowing crabs Helice crassa Dana (Grapsidae) and Macrophthalmus hirtipes (Jacquinot) (Ocypodidae). Journal of Experimental Marine Biology and Ecology, 56, 49-62.
Kamvar, Z. N., Tabima, J. F. & Grunwald, N. J. 2014. Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ, 2, e281.
Knutsen, H., Jorde, P. E., André, C. & Stenseth, N. C. 2003. Fine-scaled geographical population structuring in a highly mobile marine species: The Atlantic cod. Molecular Ecology, 12, 385-394.
Koolkalya, S., Thapanand, T., Tunkijjanujij, S., Havanont, V. & Jutagate, T. 2006. Aspects in spawning biology and migration of the mud crab Scylla olivacea in the Andaman Sea, Thailand. Fisheries Management and Ecology, 13, 391-397.
Kvingedal, R., Evans, B. S., Lind, C. E., Taylor, J. J. U., Dupont-Nivet, M. & Jerry, D. R. 2010. Population and family growth response to different rearing location, heritability estimates and genotype×environment interaction in the silver-lip pearl oyster (Pinctada maxima). Aquaculture, 304, 1-6.
Lal, M. M., Southgate, P. C., Jerry, D. R., Bosserelle, C. & Zenger, K. R. 2017. Swept away: Ocean currents and seascape features influence genetic structure across the 18,000 Km Indo-Pacific distribution of a marine invertebrate, the black-lip pearl oyster Pinctada margaritifera. BMC Genomics, 18.
Lal, M. M., Southgate, P. C., Jerry, D. R. & Zenger, K. R. 2016. Fishing for divergence in a sea of connectivity: The utility of ddRADseq genotyping in a marine invertebrate, the black-lip pearl oyster Pinctada margaritifera. Mar Genomics, 25, 57-68.
Lebata, M., Hazel, J., Vay, L. L., Primavera, J. H., Walton, M. E. & Biñas, J. B. 2007. Baseline assessment of fisheries for three species of mud crabs (Scylla spp.) in the mangroves of Ibajay, Aklan, Philippines. Bulletin of Marine Science, 80, 891-904.
Lee, S. R., Jo, Y. S., Park, C. H., Friedman, J. M. & Olson, M. S. 2018. Population genomic analysis suggests strong influence of river network on spatial distribution of genetic variation in invasive saltcedar across the southwestern United States. Mol Ecol, 27, 636-646.
Legendre, P. & Gallagher, E. D. 2001. Ecologically meaningful transformations for ordination of species data. Oecologia, 129, 271-280.
Legendre, P., Oksanen, J. & Ter Braak, C. J. F. 2011. Testing the significance of canonical axes in redundancy analysis. Methods in Ecology and Evolution, 2, 269-277.
Long, J. B. & Giri, C. 2011. Mapping the Philippines' mangrove forests using Landsat imagery. Sensors (Basel), 11, 2972-81.
Lotterhos, K. E. & Whitlock, M. C. 2014. Evaluation of demographic history and neutral parameterization on the performance of FST outlier tests. Mol Ecol, 23, 2178-92.
Lourie, S. A., Green, D. M. & Vincent, A. C. 2005. Dispersal, habitat differences, and comparative phylogeography of Southeast Asian seahorses (Syngnathidae: Hippocampus). Mol Ecol, 14, 1073-94.
Ma, H., Ma, C. & Ma, L. 2012. Molecular identification of genus Scylla (Decapoda: Portunidae) based on DNA barcoding and polymerase chain reaction. Biochemical Systematics and Ecology, 41, 41-47.
Meynecke, J.-O., Lee, S. Y., Grubert, M., Brown, I., Montgomery, S., Johnston, D. & Gillson, J. 2008. Evaluating the environmental drivers of mud crab (Scylla serrata) catches in Australia. Fisheries Research and Development Corporation and Griffith University, 1-94.
Morgulis, A., Coulouris, G., Raytselis, Y., Madden, T. L., Agarwala, R. & Schaffer, A. A. 2008. Database indexing for production MegaBLAST searches. Bioinformatics, 24, 1757-64.
Moritz, C. 1994. Defining ‘Evolutionarily Significant Units’ for conservation. Trends in Ecology & Evolution, 9, 373-375.
Moser, S., Macintosh, D., Laoprasert, S. & Tongdee, N. 2005. Population ecology of the mud crab Scylla olivacea: a study in the Ranong mangrove ecosystem, Thailand, with emphasis on juvenile recruitment and mortality. Fisheries Research, 71, 27-41.
Motoh, H., De La Peña, D. & Tampos, E. 1977. Laboratory breeding of the mud crab Scylla serrata (Forskal) through the zoea and megalopa stages to the crab stage. SEAFDEC Aquaculture Department Quarterly Research Report, 1, 14-18.
Narum, S. R. & Hess, J. E. 2011. Comparison of F(ST) outlier tests for SNP loci under selection. Mol Ecol Resour, 11 Suppl 1, 184-94.
Nayfa, M. G. & Zenger, K. R. 2016. Unravelling the effects of gene flow and selection in highly connected populations of the silver-lip pearl oyster (Pinctada maxima). Mar Genomics, 28, 99-106.
Nurdiani, R. & Zeng, C. 2007. Effects of temperature and salinity on the survival and development of mud crab, Scylla serrata (Forsskål), larvae. Aquaculture Research, 38, 1529-1538.
Obata, Y., Imai, H., Kitakado, T., Hamasaki, K. & Kitada, S. 2006. The contribution of stocked mud crabs Scylla paramamosain to commercial catches in Japan, estimated using a genetic stock identification technique. Fisheries Research, 80, 113-121.
Ogawa, C. Y., Hamasaki, K., Kitada, S., Obata, Y. & Dan, S. 2012. Species composition, reproduction, and body size of mud crabs, Scylla spp., caught in Urado Bay, Japan. Journal of Crustacean Biology, 32, 762-768.
Oksanen, J., Blanchet, F., Friendly, M., Kindt, R., Legendre, P., Mcglinn, D., Minchin, P., O’hara, R., Simpson, G. & Solymos, P. 2019. vegan: Community Ecology Package. R package version 2.5-6. 2019.
Oppo, D. W., Linsley, B. K., Rosenthal, Y., Dannenmann, S. & Beaufort, L. 2003. Orbital and suborbital climate variability in the Sulu Sea, western tropical Pacific. Geochemistry, Geophysics, Geosystems, 4, 1-20.
Palumbi, S. R. 2003. Population Genetics, Demographic Connectivity, and the Design of Marine Reserves. Ecological Applications, 13, 146-158.
Paradis, E. 2010. pegas: an R package for population genetics with an integrated-modular approach. Bioinformatics, 26, 419-20.
Paran, F. J. M. & Ravago-Gotanco, R. J. Genetic diversity and stock delineation of Philippine populations of the orange mud crab, Scylla olivacea.  Philippines: In the forefront of the mud crab industry development: proceedings of the 1st National Mud Crab Congress, 16-18 November 2015, Iloilo City, Philippines, 2017. Aquaculture Department, Southeast Asian Fisheries Development Center, 138.
Paris, C. B., Helgers, J., Van Sebille, E. & Srinivasan, A. 2013. Connectivity Modeling System: A probabilistic modeling tool for the multi-scale tracking of biotic and abiotic variability in the ocean. Environmental Modelling & Software, 42, 47-54.
Pata, P. R. & Yniguez, A. T. 2019. Larval connectivity patterns of the North Indo-West Pacific coral reefs. PLoS One, 14, e0219913.
Paterno, M., Schiavina, M., Aglieri, G., Ben Souissi, J., Boscari, E., Casagrandi, R., Chassanite, A., Chiantore, M., Congiu, L., Guarnieri, G., Kruschel, C., Macic, V., Marino, I. a. M., Papetti, C., Patarnello, T., Zane, L. & Melià, P. 2017. Population genomics meet Lagrangian simulations: Oceanographic patterns and long larval duration ensure connectivity among Paracentrotus lividus populations in the Adriatic and Ionian seas. Ecology and Evolution, 7, 2463-2479.
Peterson, B. K., Weber, J. N., Kay, E. H., Fisher, H. S. & Hoekstra, H. E. 2012. Double Digest RADseq: An Inexpensive Method for De Novo SNP Discovery and Genotyping in Model and Non-Model Species. PloS one, 7, e37135.
Quinn, N. & Kojis, B. 1987. Reproductive Biology of Scylla Spp. (Crustacea: Portunidae) from the Labu Estuary in Papua New Guinea. Bulletin of Marine Science, 41, 234-241.
Ravago-Gotanco, R. & Kim, K. M. 2019. Regional genetic structure of sandfish Holothuria (Metriatyla) scabra populations across the Philippine archipelago. Fisheries Research, 209, 143-155.
Raynal, J. M., Crandall, E. D., Barber, P. H., Mahardika, G. N., Lagman, M. C. & Carpenter, K. E. 2014. Basin isolation and oceanographic features influencing lineage divergence in the humbug damselfish <I>(Dascyllus aruanus)</I> in the Coral Triangle. Bulletin of Marine Science, 90, 513-532.
Reeves, P. A., Bowker, C. L., Fettig, C. E., Tembrock, L. R. & Richards, C. M. 2016. Effect of error and missing data on population structure inference using microsatellite data. bioRxiv, 80630-80630.
Reynolds, R. W., Smith, T. M., Liu, C., Chelton, D. B., Casey, K. S. & Schlax, M. G. 2007. Daily high-resolution-blended analyses for sea surface temperature. Journal of Climate, 20, 5473-5496.
Riginos, C., Hock, K., Matias, A. M., Mumby, P. J., Van Oppen, M. J. H. & Lukoschek, V. 2019. Asymmetric dispersal is a critical element of concordance between biophysical dispersal models and spatial genetic structure in Great Barrier Reef corals. Diversity and Distributions, 25, 1684-1696.
Riginos, C. & Liggins, L. 2013. Seascape Genetics: Populations, Individuals, and Genes Marooned and Adrift. Geography Compass, 7, 197-216.
Robertson, W. D. & Kruger, A. 1994. Size at Maturity, Mating and Spawning in the Portunid Crab Scylla serrata (Forskål) in Natal, South Africa. Estuarine, Coastal and Shelf Science, 39, 185-200.
Rochette, N. C., Rivera-Colon, A. G. & Catchen, J. M. 2019. Stacks 2: Analytical methods for paired-end sequencing improve RADseq-based population genomics. Mol Ecol, 28, 4737-4754.
Rosly, H. A., Nor, S. A., Yahya, K. & Naim, D. M. 2013. Mitochondrial DNA diversity of mud crab Scylla olivacea (Portunidae) in Peninsular Malaysia: a preliminary assessment. Mol Biol Rep, 40, 6407-18.
Sandoval-Castillo, J., Robinson, N. A., Hart, A. M., Strain, L. W. S. & Beheregaray, L. B. 2018. Seascape genomics reveals adaptive divergence in a connected and commercially important mollusc, the greenlip abalone (Haliotis laevigata), along a longitudinal environmental gradient. Mol Ecol, 27, 1603-1620.
Sanford, E. & Kelly, M. W. 2011. Local adaptation in marine invertebrates. Annual review of marine science, 3, 509-535.
Schiavina, M., Marino, I. a. M., Zane, L. & Melià, P. 2014. Matching oceanography and genetics at the basin scale. Seascape connectivity of the Mediterranean shore crab in the Adriatic Sea. Molecular Ecology, 23, 5496-5507.
Schunter, C., Carreras-Carbonell, J., Macpherson, E., Tintoré, J., Vidal-Vijande, E., Pascual, A., Guidetti, P. & Pascual, M. 2011. Matching genetics with oceanography: Directional gene flow in a Mediterranean fish species. Molecular Ecology, 20, 5167-5181.
Selkoe, K. A., D’aloia, C. C., Crandall, E. D., Iacchei, M., Liggins, L., Puritz, J. B., Von Der Heyden, S. & Toonen, R. J. 2016. A decade of seascape genetics: contributions to basic and applied marine connectivity. Marine Ecology Progress Series, 554, 1-19.
Sjöqvist, C., Godhe, A., Jonsson, P. R., Sundqvist, L. & Kremp, A. 2015. Local adaptation and oceanographic connectivity patterns explain genetic differentiation of a marine diatom across the North Sea-Baltic Sea salinity gradient. Molecular Ecology, 24, 2871-2885.
Swearer, S. E., Treml, E. A. & Shima, J. S. 2019. A Review of Biophysical Models of Marine Larval Dispersal. In: HAWKINS, S. J., ALLCOCK, A. L., BATES, A. E., FIRTH, L. B., SMITH, I. P., SWEARER, S. E. & TODD, P. A. (eds.) Oceanography and Marine Biology: An Annual Review. Taylor & Francis.
Teske, P. R., Sandoval-Castillo, J., Van Sebille, E., Waters, J. & Beheregaray, L. B. 2016. Oceanography promotes self-recruitment in a planktonic larval disperser. Scientific Reports, 6.
Thirunavukkarasu, N., Sheeba Anitha Nesacumari, C. & A. Shanmugam, A. 2014. Larval rearing and seed production of mud crab Scylla tranquebarica (Fabricius , 1798). International Journal of Fisheries and Aquatic Studies, 2, 19-25.
Tian, C., Kosoy, R., Nassir, R., Lee, A., Villoslada, P., Klareskog, L., Hammarstrom, L., Garchon, H. J., Pulver, A. E., Ransom, M., Gregersen, P. K. & Seldin, M. F. 2009. European population genetic substructure: further definition of ancestry informative markers for distinguishing among diverse European ethnic groups. Mol Med, 15, 371-83.
Toonen, R. J., Puritz, J. B., Forsman, Z. H., Whitney, J. L., Fernandez-Silva, I., Andrews, K. R. & Bird, C. E. 2013. ezRAD: a simplified method for genomic genotyping in non-model organisms. PeerJ, 1, e203.
Truelove, N. K., Kough, A. S., Behringer, D. C., Paris, C. B., Box, S. J., Preziosi, R. F. & Butler, M. J. I. 2017. Biophysical connectivity explains population genetic structure in a highly dispersive marine species. Coral Reefs, 36, 233-244.
Van Wyngaarden, M., Snelgrove, P. V., Dibacco, C., Hamilton, L. C., Rodriguez-Ezpeleta, N., Jeffery, N. W., Stanley, R. R. & Bradbury, I. R. 2017. Identifying patterns of dispersal, connectivity and selection in the sea scallop, Placopecten magellanicus, using RADseq-derived SNPs. Evol Appl, 10, 102-117.
Van Wyngaarden, M., Snelgrove, P. V. R., Dibacco, C., Hamilton, L. C., Rodriguez-Ezpeleta, N., Zhan, L., Beiko, R. G. & Bradbury, I. R. 2018. Oceanographic variation influences spatial genomic structure in the sea scallop, Placopecten magellanicus. Ecol Evol, 8, 2824-2841.
Vilas, A., Perez-Figueroa, A. & Caballero, A. 2012. A simulation study on the performance of differentiation-based methods to detect selected loci using linked neutral markers. J Evol Biol, 25, 1364-76.
Viswanathan, C., Pravinkumar, M., Suresh, T. V., Elumalai, V. & Raffi, S. M. 2019. Reproductive biology of the orange mud crab Scylla olivacea (Herbst, 1796) from the Pichavaram mangroves of south-east India. Indian Journal of Fisheries, 66.
Wang, L., Liu, S., Zhuang, Z., Guo, L., Meng, Z. & Lin, H. 2013. Population genetic studies revealed local adaptation in a high gene-flow marine fish, the small yellow croaker (Larimichthys polyactis). PLoS One, 8, e83493.
Wang, S., Meyer, E., Mckay, J. K. & Matz, M. V. 2012. 2b-RAD: a simple and flexible method for genome-wide genotyping. Nat Methods, 9, 808-10.
Waples, R. S. 2010. Spatial-temporal stratifications in natural populations and how they affect understanding and estimation of effective population size. Mol Ecol Resour, 10, 785-96.
Warnes, G., Bolker, B., Bonebakker, L., Gentleman, R., Liaw, W., Lumley, T., Maechler, M., Magnusson, A., Moeller, S. & Schwartz, M. 2020. gplots: Various R Programming Tools for Plotting Data. R package version 3.1.0.
Warnes, G., Gorjanc, G., Leisch, F. & Man, M. 2019. genetics: Population genetics. R package version 1.3. 8.1. 2. .
Weir, B. S. & Cockerham, C. C. 1984. Estimating F-Statistics for the Analysis of Population Structure. Evolution, 38, 1358-1370.
Whitlock, M. C. & Mccauley, D. E. 1999. Indirect measures of gene flow and migration: FST not equal to 1/(4Nm + 1). Heredity (Edinb), 82 ( Pt 2), 117-25.
Xuereb, A., Benestan, L., Normandeau, E., Daigle, R. M., Curtis, J. M. R., Bernatchez, L. & Fortin, M. J. 2018. Asymmetric oceanographic processes mediate connectivity and population genetic structure, as revealed by RADseq, in a highly dispersive marine invertebrate (Parastichopus californicus). Mol Ecol, 27, 2347-2364.
Zhang, Z., Schwartz, S., Wagner, L. & Miller, W. 2000. A greedy algorithm for aligning DNA sequences. J Comput Biol, 7, 203-14.
 

Tables (separate file)

Table 1. Summary information of S. olivacea sample location and genetic diversity estimates.
 
Table 2. Filtering processes and specific parameters used to filter the final SNP panel of S. olivacea. Each filtering step indicates the number of SNPs that were removed and retained. Filtering tools and programs are indicated in the methods section.
 
Table 3. Summary of the redundancy analysis (RDA) of genetic data (neutral and outlier loci-only) and environmental variables including the asymmetric eigenvector maps (AEMs), sea surface temperature (SST), and rainfall. Only those variables retained by backward selection and were significant (p < 0.05) are included in the final model. R2adj 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.
 

Figure legends

Figure 1. Estimated larval dispersal and connectivity of S. olivacea populations in the Sulu Sea basin from Lagrangian simulations. (a) Node connections plot generated from the larval dispersal estimates resulted in six population clusters indicated by the colored nodes and connections, gray links indicating connections between clusters. Site codes represent particles release sites and settlement nodes were labelled numerically. (b) Population connectivity matrix showing the proportion of larvae successfully settled from the source locations (y-axis) to settlement areas (x-axis). Sulu Sea populations were segregated according to their boundary location (east or west). See Table 1 for location codes.
 
Figure 2. Spatial patterns of S. olivacea connectivity based on (a) empirical genetic estimates (pairwise FST) from the putatively neutral loci (1,643 SNPs); and (b) asymmetric larval dispersal estimates (AEM6) from the particle dispersal simulations in the Sulu Sea. Sampling points are colored according to their cluster assignment. Persistent surface currents in the Sulu Sea are shown.
 
Figure 3. Heatmap of pairwise genetic differentiation (FST) of S. olivacea between sampling locations in the Sulu Sea, based on Weir and Cockerham weighted estimates using (a) all markers (1,655 loci), (b) neutral-only markers (1,643 SNPs), and (c) outlier only (12 SNPs). All points were clustered by pairwise FST values, based on classification and hierarchical clustering method. Site combinations (below diagonal) with boxed bold lines indicate significant genetic differentiation following FDR correction for multiple tests (p < 0.05).
 
Figure 4. (a) Scatterplot of genetic clusters of Sulu Sea populations identified by DAPC based on neutral loci (1,643 SNPs) and sampling location information as a prior, with density plots for (b) the first and (c) the second discriminant axes. Genetic clusters identified using outlier loci (12 SNPs) and sampling location information as a prior (d).
 
Figure 5. Distribution of genetic clusters of S. olivacea in the Sulu Sea based on outlier loci and Geneland analysis. 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.
 

Data Accessibility Statement

Raw demultiplexed sequence libraries in fastq format were archived in NCBI SRA (BioProject Accession # PJRNA662443) and will be released upon publication. The corresponding filtered datasets and R scripts for analyses have been deposited in Dryad Digital Repository: https://doi.org/10.5061/dryad.3xsj3txdz
 

Author Contribution Section

MJM: Conceptualization (equal), Data Curation, Formal Analysis (lead), Investigation, Methodology (equal), Visualization, Writing-Original Draft (equal)
RRG: Conceptualization (equal), Formal Analysis (supporting), Funding Acquisition, Methodology (equal), Project Administration, Supervision, Writing-Original Draft (equal)
 

Acknowledgements

This project was funded by the Department of Science and Technology - Philippine Council for Agriculture, Aquatic, and Natural Resources Research and Development (DOST-PCAARRD) and implemented by the University of the Philippines Marine Science Institute (UP-MSI). We are deeply thankful to Dr. Evangeline Magdaong, Jeniffer De Maligaya, and Benedict Castro of the Physical Oceanography Laboratory, UP-MSI headed by Dr. Cesar Villanoy for the larval dispersal biophysical modelling. We are grateful to Angela Camille Aguila and Simon Alcantara for laboratory assistance, Bhenjamin Ona for the remote sensing data, Dr. Din Matias for analysis recommendations, Von Yip for QGIS assistance, Dr. Richard Mualil and Yunadzmal Ong of Mindanao State University (MSU Tawi-Tawi) for sample collection. We thank Sharon Magnuson and Chris Bird (Genomics Core Lab, Texas A&M University, Corpus Christi) for RAD sequencing.
 

Ethical Approval

Scylla olivacea is a commercially harvested species and was not considered as a regulated species during the time of sampling. For Palawan, collection and local transport were covered under Palawan Council for Sustainable Development (PCSD) local transport permits and Wildlife Gratuitous Permit (GP) no. 2016-23. Export of material for DNA sequencing was covered under the Bureau of Fisheries and Aquatic Resources GP no. 2018-0005.
 

Conflict of Interest

None declared.