\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,english]{babel}



\usepackage{float}








\begin{document}

\title{Genetic differentiation and signatures of local adaptation revealed by
RADseq for a highly-dispersive mud crab~\emph{Scylla olivacea} (Herbst,
1796) in the Sulu Sea}



\author[1]{Michael John  Mendiola}%
\author[1]{Rachel Ravago-Gotanco}%
\affil[1]{The Marine Science Institute, University of the Philippines Diliman, Quezon City 1101 Philippines}%


\vspace{-1em}



  
  \date{May 26, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Connectivity of marine populations is shaped by complex interactions of
biological and physical processes across the seascape. The influence of
environmental features on the genetic structure of populations has key
implications to 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 genetics
approach combining larval dispersal modeling and population genomic
analysis based on RADseq to examine environmental factors influencing
patterns of genetic structure and connectivity for a highly-dispersive
mud crab, \emph{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,
significant genetic differentiation was detected among eight Sulu Sea
populations based on 1,655 single-nucleotide polymorphism (SNP) markers
(\emph{F}\textsubscript{ST} = 0.0057, \emph{p} = 0.001) and a subset of
1,643 putatively neutral SNP markers (\emph{F}\textsubscript{ST} =
0.0042, \emph{p} = 0.001). Oceanography influences genetic structure, as
redundancy analysis (RDA) revealed significant contribution of
asymmetric ocean currents to neutral genetic variation
(\emph{R}\textsuperscript{2}\textsubscript{adj} = 0.133; \emph{p} =
0.035). Genetic structure may also reflect demographic factors, with
divergent populations characterized by low effective population sizes
(\emph{N}\textsubscript{e} \textless{} 50). Pronounced latitudinal
genetic structure was recovered for loci putatively under selection
(\emph{F}\textsubscript{ST} = 0.2390, \emph{p} = 0.001), significantly
correlated with variability in mean sea surface temperatures during peak
spawning months of \emph{S. olivacea}
(\emph{R}\textsuperscript{2}\textsubscript{adj} = 0.763; \emph{p} =
0.041), suggesting putative signatures of selection and local adaptation
of early life history stages to thermal clines. 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.%
\end{abstract}%



\sloppy


\par\null

\section*{1 Introduction}\label{introduction}

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 and Sponaugle, 2009). The
absence of apparent physical barriers in the ocean combined with high
dispersal potentials characteristic of most marine organisms shaped the
paradigm of genetic homogeneity in the marine environment (Hauser and
Carvalho, 2008). Advances in DNA sequencing technologies now provide the
ability to interrogate genetic variation genome-wide, providing adequate
resolution to study population genetic processes (Davey et al., 2011,
Andrews et al., 2016). High-throughput genotyping of single nucleotide
polymorphisms (SNPs) employing restriction-site associated DNA
sequencing approaches (RAD-sequencing; Baird et al., 2008) is a
widely-used approach having several variants, e.g ddRAD (Peterson et
al., 2012), 2b-RAD (Wang et al., 2012), and ezRAD (Toonen et al., 2013).
Genetic approaches, coupled with increased capabilities in numerical
modeling to simulate and track larval dispersal in the marine
environment (reviewed in Swearer et al., 2019), have resulted in
considerable interest and wide use of seascape genomics approaches to
examine the processes that shape genetic variation in the marine
environment (Riginos and Liggins, 2013, 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) allows the
examination of environmental factors and their influence on genetic
structure. There is a growing body of literature documenting genetic
structure either due to neutral variation or adaptive polymorphisms, at
finer spatial scales than expected from species dispersal potentials.
Recent studies have reported environmental factors influencing neutral
and adaptive genetic variation, among them: ocean currents (Benestan et
al., 2016, Gilg and 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 ~(Wang et al. 2013; Chu et al. 2014,
Sandoval-Castillo et al. 2018, Hoey and Pinsky 2018), and salinity
(Sj\selectlanguage{ngerman}öqvist et al., 2015).

Mud crabs (genus: \emph{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 occur in the
Philippines (\emph{S. olivacea, S. tranquebarica}, and \emph{S.
serrata}), with \emph{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) of 20 to 30
days (Jantrarotai et al., 2002, Motoh et al., 1977, Thirunavukkarasu et
al., 2014), which can extend up to 75 days depending on environmental
conditions (Baylon, 2010). Although ocean currents may play a huge role
in larval dispersal and settlement success (Cowen and Sponaugle, 2009),
survival and development of mud crab larvae strongly depends on water
temperature and salinity (Baylon, 2011, Hamasaki, 2003, Hill, 1974,
Nurdiani and Zeng, 2007). In the Sulu Sea basin, \emph{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 westward-flowing 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 \emph{S. olivacea} such as the seahorse \emph{Hippocampus
spinosissimus} (Lourie et al., 2005), damselfish \emph{Dascyllus
aruanus} (Raynal et al., 2014) and sea cucumber \emph{Holuthuria scabra}
(Ravago-Gotanco and Kim, 2019). There is limited information however, on
the genetic structure of Philippine populations of \emph{S. olivacea,}
with one study reporting weak but significant genetic differentiation of
Philippine populations based on microsatellite loci (Paran and
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 (\emph{S. olivacea}) in the Sulu Sea basin. Using a
seascape genomics approach, we aimed to (1) characterize genetic
structure of \emph{S. olivacea} across the Sulu Sea and examine spatial
patterns of genetic connectivity using SNP markers generated by RADseq;
(2) examine the influence of oceanographic circulation on genetic
structure and connectivity of \emph{S. olivacea} populations in the Sulu
Sea; and (3) 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 \emph{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 rainfall
(as a proxy for salinity) on adaptive divergence of \emph{S. olivacea}.
This study provides valuable insights into the spatial scales of
dispersal, patterns of genetic structure, and the influence of
environmental and evolutionary processes on population connectivity of
\emph{S. olivacea} in the Sulu Sea basin. The results provide
information useful to support the development of management and
conservation strategies for the fishery resource.

\par\null

\section*{2 Materials and methods}\label{materials-and-methods}

\subsection*{2.1 Larval dispersal simulation in the Sulu Sea
basin}\label{larval-dispersal-simulation-in-the-sulu-sea-basin}

A larval dispersal model was developed to examine connectivity of
\emph{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°
(\textasciitilde{}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 \emph{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 \emph{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
\emph{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 \emph{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
(\textasciitilde{}75 km radius), resulting in a reduced population
source-sink dataset.

~

\subsection*{2.2 Sample collection, species identification, and total DNA
extraction}\label{sample-collection-species-identification-and-total-dna-extraction}

Adult \emph{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\selectlanguage{english} spectrophotometer.

~

\subsection*{2.3 Double-digest RAD (ddRAD)
sequencing}\label{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 \emph{MluCI} and \emph{MspI} and \emph{}
purified using AMPureXP (Beckman Coulter). Custom barcoded adapters P1
and P2 (see Peterson et al. 2012 for sequences) were then ligated to
\textasciitilde{}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\selectlanguage{english} 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.

~

\subsection*{2.4 Read processing and SNP
filtering}\label{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
\emph{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 \emph{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 \emph{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 \emph{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 \emph{pegas} v0.11 (Paradis, 2010),
excluding loci which exhibited significant deviation from HWE (p
\textless{} 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 \emph{genetics} v1.3.8.1.2 (Warnes et al., 2019), and SNPs
at strong linkage disequilibrium (\emph{r}\textsuperscript{2}
\textgreater{} 0.8) were removed (Lee et al., 2018, Tian et al., 2009).
SNPs with high observed heterozygosity (H\textsubscript{o}
\textgreater{} 0.6) were also dropped from the dataset (e.g. (Ackiss et
al., 2018, Hohenlohe et al., 2010, Van Wyngaarden et al., 2017) using
\emph{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).

~

\subsection*{2.5 Identifying non-neutral
SNPs}\label{identifying-non-neutral-snps}

Loci were identified as being putatively neutral or under selection
using two differentiation-based (\emph{F}\textsubscript{\emph{ST}})
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
\emph{F}\textsubscript{ST} to the observed global
\emph{F}\textsubscript{ST} 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) \emph{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 \emph{F}\textsubscript{ST}
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).

~

\subsection*{2.6 Population genetic structure and effective population
size}\label{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
\emph{F}\textsubscript{ST} (1984) to estimate genetic differentiation
over all populations (global \emph{F}\textsubscript{ST}) and between
populations (pairwise \emph{F}\textsubscript{ST}), calculated using the
R packages \emph{hierfstat} v0.04-22 (Goudet and Jombart, 2019) and
\emph{dartR} v1.1.11 (Gruber et al., 2018) respectively. Significance of
\emph{F}\textsubscript{ST} values was tested with 10,000 bootstrap
replicates, and p\selectlanguage{english}-values for population pairwise comparisons were
adjusted for multiple tests using the false discovery rate (FDR).
Pairwise \emph{F}\textsubscript{ST} values were visualized using
heatmaps with dendrograms generated from hierarchical clustering
analysis performed using the function \emph{heatmap.2} from the
\emph{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 \emph{adegenet} v2.1.3 (Jombart, 2008).
DAPC was performed using sampling location as prior information.
Cross-validation was performed using the function \emph{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\selectlanguage{english}-in of 200, maximum rate of the Poisson process fixed to 146 (N
= number of individuals), maximum number of nuclei in the
Poisson\selectlanguage{english}-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
\emph{F}\textsubscript{ST} analysis and DAPC on the three SNP datasets.
AMOVA was performed using \emph{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 \emph{gl.ibd} function in \emph{dartR} which performs a Mantel test
(number of permutations = 9,999) to assess correlation of log(geographic
distance) and linearized genetic distance
(\emph{F}\textsubscript{ST}/(1-\emph{F}\textsubscript{ST})) matrices.
Geographic distances were measured as the shortest distance over water
between all pairs of sites in the Sulu Sea using the \emph{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
(\emph{N}\textsubscript{\emph{eLD}}) of NeEstimator v2.01 ~(Do et al.,
2014).

~

\subsection*{2.7 Genetic structure and environmental
factors}\label{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
\emph{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 \emph{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.25\textsuperscript{o} x
0.25\textsuperscript{o}) monthly mean data in the Sulu Sea domain
(5-14\textsuperscript{o}N, 116-124\textsuperscript{o}E) from 1987-2005
were downloaded from the NOAA/OAR/ESRL PSL website at
\url{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
\emph{ordistep} in \emph{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
(\emph{R}\textsuperscript{2}\textsubscript{adj}) was calculated. Partial
RDA was performed, and the global analysis of variance was determined
using the \emph{anova} function in \emph{stats v3.5.3} with 999
permutations. Significance of individual RDA axes were also assessed
using \emph{anova} (permutations = 999) and selected environmental
vectors were fitted into the ordination using \emph{envfit}
(permutations = 999).

~

\section*{3 Results}\label{results}

\subsection*{3.1 Larval dispersal
estimates}\label{larval-dispersal-estimates}

Lagrangian simulations of larval dispersal using a parameterized
biophysical model demonstrate ocean current-mediated connectivity among
\emph{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.

~

\subsection*{3.2 SNP filtering and identification of non-neutral
loci}\label{snp-filtering-and-identification-of-non-neutral-loci}

~~~~~~~~~~~ A total of 661,372,129 paired-end (PE) reads from 146
individual \emph{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
\emph{F}\textsubscript{ST} 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
\emph{F}\textsubscript{ST} \emph{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 (\emph{Portunus trituberculatus, Penaeus vannamei})
and fish (\emph{Chanos chanos, Danio renio}). One tag had 87.76\%
identity match to a microsatellite region of Japanese blue swimming crab
\emph{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.

~

\subsection*{~3.3 Genetic differentiation and population
structure}\label{genetic-differentiation-and-population-structure}

All~\emph{Scylla olivacea} populations exhibited significant genetic
differentiation based on global estimates of~\emph{F}\textsubscript{ST}
calculated using all loci (1655 SNPs;~\emph{F}\textsubscript{ST} =
0.0070,~\emph{p} = 0.001) and putatively neutral loci (1643
SNPs;~\emph{F}\textsubscript{ST} = 0.0056, p = 0.001). Excluding
outgroup sites, Sulu Sea samples (n = 116 individuals; 8 sites) still
exhibited significant genetic differentiation for all loci
(\emph{F}\textsubscript{ST} = 0.0057,~\emph{p} = 0.001) and putatively
neutral loci (\emph{F}\textsubscript{ST} = 0.0042,~\emph{p} = 0.001).~
Pairwise~\emph{F}\textsubscript{ST} values using all loci revealed the
most genetically divergent populations to be CRN
(\emph{F}\textsubscript{ST} range = 0.0021 - 0.0219), PPC
(\emph{F}\textsubscript{ST} range = 0.0014 - 0.0187), and one outgroup
GSC (\emph{F}\textsubscript{ST} range = 0.0014 - 0.0219) (Table S3). The
divergence of CRN, PPC and GSC was also evident for neutral loci,
although~\emph{F}\textsubscript{\emph{ST}} estimates were slightly lower
(Table S4). Dendrograms of pairwise~\emph{F}\textsubscript{ST} values
clearly separate CRN and PPC from the rest of the Sulu Sea populations
using all loci and putatively neutral loci; all~\emph{p}-values for
pairwise comparisons were \textless{} 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
(\emph{F}\textsubscript{ST} = 0.0016,~\emph{p} = 0.002), with
significant~\emph{F}\textsubscript{ST} 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 \emph{F}\textsubscript{ST} \emph{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 (\emph{F}\textsubscript{\emph{CT}} =
0.017, \emph{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, \emph{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, \emph{p} = 0.083).
Using the putatively neutral loci dataset, genetically divergent
populations CRN, PPC and GSC were estimated to have small effective
population size (\emph{N}\textsubscript{e}: CRN = 24.6-26.8, PPC =
10.7-11.2, GSC = 9.6-10.1) compared to other localities where
\emph{N}\textsubscript{e} ~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; \emph{F}\textsubscript{ST} = 0.2390, \emph{p} = 0.001). A
dendrogram based on pairwise \emph{F}\textsubscript{ST} 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
\emph{F}\textsubscript{ST} among populations within each of the 3
clusters are not significant (FDR adjusted \emph{p}~ \textgreater{}
0.05), i.e. for ROX-PPC-ANT and BAT-TWI, while between-cluster
comparisons are significant (FDR adjusted \emph{p} \textless{} 0.05)
(Table S5). \textsubscript{} 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 \emph{F}\textsubscript{ST}, 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
(\emph{F}\textsubscript{CT} = 0.209; \emph{p} = 0.001) accounting for
20.9\% of the total observed variance.~ No further structure is detected
among samples within groups (\emph{F}\textsubscript{SC} = 0.0168 ,
\emph{p} = 0.166).~ However, the discordance between pairwise
\emph{F}\textsubscript{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
\emph{F}\textsubscript{ST} and multivariate methods. To examine this
further, we performed two separate analyses to handle missing data: (1)
remove genotypes (individuals) with \textgreater{} 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 12\textsuperscript{o}N),
ANT-PPC-ROX (at 9\textsuperscript{o}N - 11\textsuperscript{o}N), as
genetically distinct groups from BAT-NEG and TWI (at
5\textsuperscript{o}N to 8\textsuperscript{o}N) by pairwise
\emph{F}\textsubscript{\emph{ST}}, DAPC and GENELAND, provides support
for a pattern of latitudinal structure of Sulu Sea populations based on
outlier loci.

~

\subsection*{3.4 Genetic structure and environmental
factors}\label{genetic-structure-and-environmental-factors-1}

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
\emph{p} \textgreater{} 0.05) (Table 4). Together, these two AEM
eigenfunctions explained 13.3\% of neutral genetic variation among sites
(adjusted \emph{R}\textsuperscript{2}; \emph{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
(\emph{p} = 0.065), whereas RDA2 accounted for 39.1\% of the total
genetic variation (\emph{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 (\emph{p}
\textless{} 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 \emph{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
(\emph{R}\textsuperscript{2}\textsubscript{adj} \textsuperscript{} =
0.763; \emph{p} = 0.041). RDA1 explained the highest fraction of genetic
variation comprising 88.0\% (\emph{p} = 0.048), while RDA2 accounted for
10.2\% (\emph{p} = 0.341). Moreover, rainfall data explained a lower
proportion of the variation
(\emph{R}\textsuperscript{2}\textsubscript{adj} = 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 (\emph{p} = 0.089), which suggests that mean SST
(Figure S2) is a stronger predictor of the observed latitudinal genetic
variation than rainfall.

\section*{4 Discussion}\label{discussion}

In this study, we employed a seascape genomics approach to examine
environmental factors influencing genetic structure of \emph{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 \emph{S. olivacea} populations in the Sulu Sea, which may be used to
devise spatially explicit management and conservation interventions.

~

\subsection*{4.1 Genetic structure and connectivity in the Sulu
Sea}\label{genetic-structure-and-connectivity-in-the-sulu-sea}

Broad dispersal of \emph{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 \emph{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 \emph{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 \emph{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
\emph{F}\textsubscript{ST} values over sampling scales of hundreds of
kilometers have been reported for broadly-dispersing species such as
highly mobile fish (Atlantic cod, \emph{F}\textsubscript{\emph{ST}} =
0.004; Knutsen et al., 2003), or invertebrates with extensive larval
duration periods such as red rock lobsters
(\emph{F}\textsubscript{\emph{ST}} \textsubscript{} = 0.004; Iacchei et
al., 2013) and spiny lobsters (\emph{F}\textsubscript{\emph{ST}} =
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 \emph{Engraulis
encrasicolus} (Bembo et al., 1996), shore crab \emph{Carcinus aestuarii}
(Schiavina et al., 2014), and peacock wrasse \emph{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
(\emph{R}\textsuperscript{2}\textsubscript{adj} \textsuperscript{} =
0.133, \emph{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. \emph{F}\textsubscript{ST}-based approaches
reveal PPC as divergent from other Sulu Sea populations (Figure 2a, 3a,
3b). The separation of TWI, while not supported by pairwise
\emph{F}\textsubscript{ST} \emph{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 (\emph{F}\textsubscript{ST} =
0.0004, \emph{p} = 0.094), suggesting no apparent barriers to \emph{S.
olivacea} larval dispersal across the Sulu Sea. Genetic studies for
other species with similar limited adult movement, but shorter pelagic
larval durations than \emph{S. olivacea} indicate limited gene flow
between eastern and western boundary populations. Genetic structure for
the seahorse \emph{Hippocampus spinosissimus} (Lourie et al., 2005),
damselfish \emph{Dascyllus aruanus} (Raynal et al., 2014), and sea
cucumber \emph{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 (\emph{Acropora millepora}), sea urchin
(\emph{Tripneustes gratilla}), and a reef fish (\emph{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 \emph{N}\textsubscript{e} =
24.6-26.8, PPC \emph{N}\textsubscript{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,
e.g.~ red cusk-eel \emph{Genypterus chilensis} due to high fishing
pressure (Cordova-Alarcon et al., 2019), and population bottlenecks for
\emph{Gadus} morhua (Andreev et al., 2015) and \emph{Epinephelus
marginatus} (Buchholz-Sorensen 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.

~

\subsection*{4.2 Latitudinal adaptation of S. olivacea in the Sulu
Sea}\label{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). \emph{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 (\emph{F}\textsubscript{ST}, DAPC,
GENELAND), with AMOVA indicating significant differentiation among
groups accounting for 20.9\% of the total variance
(\emph{F}\textsubscript{CT} = 0.209; \emph{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 \emph{S. olivacea,} with SST
(\emph{R}\textsuperscript{2}\textsubscript{adj} \textsuperscript{} =
0.763, \emph{p} = 0.041) as a stronger predictor of genetic variation
than rainfall (\emph{R}\textsuperscript{2}\textsubscript{adj}
\textsuperscript{} = 0.656, \emph{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
\emph{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 \emph{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,
\emph{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 \emph{S. olivacea}, \emph{}
which occurs from July to November at latitudes between
9\textsuperscript{o}N to 11\textsuperscript{o}N (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 (\emph{Nucella lapillus}) along
midcoastal Maine (between 43degN and 44degN; with water and air difference
reaching up to 5-10 degC), 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 (\emph{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-2degC
thermal difference between sites (Wang et al., 2013). For a marine
diatom \emph{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 (Sjoqvist 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.

\subsection*{4.3 Implications to
Management}\label{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 \emph{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 \emph{S. olivacea}.

From the perspective of neutral loci, \emph{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 \emph{N}\textsubscript{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 (\emph{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 \emph{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 \emph{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 \emph{S. olivacea} in the
Philippines. \emph{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 \emph{S. serrata},
\emph{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.

\section*{Literature cited}\label{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.
\emph{Ecol Evol,} 8\textbf{,} 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.
\emph{Hydrobiologia,} 763\textbf{,} 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. \emph{Heliyon,} 6\textbf{,} 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). \emph{Marine Biology,} 162\textbf{,} 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. \emph{Nature Reviews Genetics,} 17\textbf{,} 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.
\emph{Am J Hum Genet,} 71\textbf{,} 304-11.

Arriola, F. J. 1940. \emph{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.
\emph{PloS one,} 3\textbf{,} 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). \emph{Journal of the World
Aquaculture Society,} 41\textbf{,} 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. \emph{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. \emph{Marine Biology,} 126\textbf{,} 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). \emph{Molecular
ecology,} 25\textbf{,} 5073-5092.

Benjamini, Y. \& Hochberg, Y. 1995. Controlling the False Discovery
Rate: A Practical and Powerful Approach to Multiple Testing.
\emph{Journal of the Royal Statistical Society. Series B
(Methodological),} 57\textbf{,} 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.). \emph{In:} RESOURCES, B.
O. F. A. A. (ed.) \emph{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.
\emph{Journal of Biogeography,} 35\textbf{,} 1593-1607.

Buchholz-Sorensen, 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. \emph{PLOS ONE,} 11\textbf{,} 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. \emph{World Congress on Aquatic Protected Areas
Proceedings}\textbf{,} 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. \emph{Sci Rep,} 7\textbf{,} 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.
\emph{Journal of Marine Systems,} 65\textbf{,} 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. \emph{Mol Ecol,} 23\textbf{,} 1863-73.

Cordova-Alarcon, 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. \emph{Frontiers in Marine Science,} 6\textbf{,} 748-748.

Cowen, R. K. \& Sponaugle, S. 2009. Larval dispersal and marine
population connectivity. \emph{Ann Rev Mar Sci,} 1\textbf{,} 443-66.

Csardi, G. \& Nepusz, T. 2006. The igraph software package for complex
network research. \emph{InterJournal, complex systems,} 1695\textbf{,}
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. \emph{Nature Reviews
Genetics,} 12\textbf{,} 499-510.

Dawson, M. N., Raskoff, K. A. \& Jacobs, D. K. 1998. Field preservation
of marine invertebrate tissue for DNA analyses. \emph{Molecular marine
biology and biotechnology,} 7\textbf{,} 145-152.

Dewoody, Y. D. \& Dewoody, J. A. 2005. On the Estimation of Genome-wide
Heterozygosity Using Molecular Markers. \emph{Journal of Heredity,}
96\textbf{,} 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. \emph{Mol Ecol Resour,} 14\textbf{,} 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: \url{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. \emph{Mol Ecol,} 16\textbf{,} 463-75.

Edmands, S. \& Timmerman, C. C. 2003. Modeling Factors Affecting the
Severity of Outbreeding Depression. \emph{Conservation Biology,}
17\textbf{,} 883-892.

Ewels, P., Magnusson, M., Lundin, S. \& Kaller, M. 2016. MultiQC:
summarize analysis results for multiple tools and samples in a single
report. \emph{Bioinformatics,} 32\textbf{,} 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. \emph{Mol Ecol Resour,} 10\textbf{,} 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.
\emph{Genetics,} 131\textbf{,} 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. \emph{Sci Rep,} 5\textbf{,} 13758.

Foll, M. \& Gaggiotti, O. 2008. A Genome-Scan Method to Identify
Selected Loci Appropriate for Both Dominant and Codominant Markers: A
Bayesian Perspective. \emph{Genetics,} 180\textbf{,} 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. \emph{Evolutionary Applications,}
8\textbf{,} 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. \emph{Aquaculture,} 173\textbf{,} 117-129.

Gilg, M. R. \& Hilbish, T. J. 2003. The Geography of Marine Larval
Dispersal: Coupling Genetics with Fine-Scale Physical Oceanography.
\emph{Ecology,} 84\textbf{,} 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. \emph{Mol Ecol Resour,} 18\textbf{,}
691-699.

Guillot, G., Mortier, F. \& Estoup, A. 2005. Geneland: a computer
package for landscape genetics. \emph{Molecular Ecology Notes,}
5\textbf{,} 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 (Forskal) (Brachyura: Portunidae) reared in the laboratory.
\emph{Aquaculture,} 219\textbf{,} 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. \emph{Dynamics of Atmospheres and Oceans,}
47\textbf{,} 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. \emph{Conserv Biol,} 25\textbf{,} 438-49.

Hauser, L. \& Carvalho, G. R. 2008. Paradigm shifts in marine fisheries
genetics: ugly hypotheses slain by beautiful facts. \emph{Fish and
Fisheries,} 9\textbf{,} 333-362.

Hill, B. J. 1974. Salinity and temperature tolerance of zoeae of the
portunid crab Scylla serrata. \emph{Marine Biology,} 25\textbf{,} 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. \emph{Am Nat,} 188\textbf{,}
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. \emph{PLoS Genet,}
6\textbf{,} 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. \emph{Journal of
hydrometeorology,} 8\textbf{,} 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/12deg and 1/25deg Global HYCOM and EAS NCOM.
\emph{Oceanography,} 24\textbf{,} 28-47.

Hyland, S., Hill, B. \& Lee, C. 1984. Movement within and between
different habitats by the portunid crab Scylla serrata. \emph{Marine
Biology,} 80\textbf{,} 57-61.

Iacchei, M., Ben-Horin, T., Selkoe, K. A., Bird, C. E.,
Garcia-Rodriguez, F. J. \& Toonen, R. J. 2013. Combined analyses of
kinship and FST suggest potential drivers of chaotic genetic patchiness
in high gene-flow populations. \emph{Molecular Ecology,} 22\textbf{,}
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. \emph{Kasetsart
Journal (Natural Science),} 36\textbf{,} 278-284.

Jombart, T. 2008. adegenet: a R package for the multivariate analysis of
genetic markers. \emph{Bioinformatics,} 24\textbf{,} 1403-5.

Jombart, T., Devillard, S. \& Balloux, F. 2010. Discriminant analysis of
principal components: a new method for the analysis of genetically
structured populations. \emph{BMC Genet,} 11\textbf{,} 94.

Jones, M. B. \& Simons, M. J. 1981. Habitat preferences of two estuarine
burrowing crabs Helice crassa Dana (Grapsidae) and Macrophthalmus
hirtipes (Jacquinot) (Ocypodidae). \emph{Journal of Experimental Marine
Biology and Ecology,} 56\textbf{,} 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. \emph{PeerJ,} 2\textbf{,} e281.

Knutsen, H., Jorde, P. E., Andre, C. \& Stenseth, N. C. 2003.
Fine-scaled geographical population structuring in a highly mobile
marine species: The Atlantic cod. \emph{Molecular Ecology,} 12\textbf{,}
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. \emph{Fisheries
Management and Ecology,} 13\textbf{,} 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
genotypexenvironment interaction in the silver-lip pearl oyster
(Pinctada maxima). \emph{Aquaculture,} 304\textbf{,} 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.
\emph{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. \emph{Mar Genomics,} 25\textbf{,} 57-68.

Lebata, M., Hazel, J., Vay, L. L., Primavera, J. H., Walton, M. E. \&
Binas, J. B. 2007. Baseline assessment of fisheries for three species of
mud crabs (Scylla spp.) in the mangroves of Ibajay, Aklan, Philippines.
\emph{Bulletin of Marine Science,} 80\textbf{,} 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. \emph{Mol Ecol,}
27\textbf{,} 636-646.

Legendre, P. \& Gallagher, E. D. 2001. Ecologically meaningful
transformations for ordination of species data. \emph{Oecologia,}
129\textbf{,} 271-280.

Legendre, P., Oksanen, J. \& Ter Braak, C. J. F. 2011. Testing the
significance of canonical axes in redundancy analysis. \emph{Methods in
Ecology and Evolution,} 2\textbf{,} 269-277.

Long, J. B. \& Giri, C. 2011. Mapping the Philippines' mangrove forests
using Landsat imagery. \emph{Sensors (Basel),} 11\textbf{,} 2972-81.

Lotterhos, K. E. \& Whitlock, M. C. 2014. Evaluation of demographic
history and neutral parameterization on the performance of FST outlier
tests. \emph{Mol Ecol,} 23\textbf{,} 2178-92.

Lourie, S. A., Green, D. M. \& Vincent, A. C. 2005. Dispersal, habitat
differences, and comparative phylogeography of Southeast Asian seahorses
(Syngnathidae: Hippocampus). \emph{Mol Ecol,} 14\textbf{,} 1073-94.

Ma, H., Ma, C. \& Ma, L. 2012. Molecular identification of genus Scylla
(Decapoda: Portunidae) based on DNA barcoding and polymerase chain
reaction. \emph{Biochemical Systematics and Ecology,} 41\textbf{,}
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. \emph{Fisheries
Research and Development Corporation and Griffith University}\textbf{,}
1-94.

Morgulis, A., Coulouris, G., Raytselis, Y., Madden, T. L., Agarwala, R.
\& Schaffer, A. A. 2008. Database indexing for production MegaBLAST
searches. \emph{Bioinformatics,} 24\textbf{,} 1757-64.

Moritz, C. 1994. Defining `Evolutionarily Significant Units' for
conservation. \emph{Trends in Ecology \& Evolution,} 9\textbf{,}
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. \emph{Fisheries Research,} 71\textbf{,} 27-41.

Motoh, H., De La Pena, D. \& Tampos, E. 1977. Laboratory breeding of the
mud crab Scylla serrata (Forskal) through the zoea and megalopa stages
to the crab stage. \emph{SEAFDEC Aquaculture Department Quarterly
Research Report,} 1\textbf{,} 14-18.

Narum, S. R. \& Hess, J. E. 2011. Comparison of F(ST) outlier tests for
SNP loci under selection. \emph{Mol Ecol Resour,} 11 Suppl 1\textbf{,}
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). \emph{Mar Genomics,} 28\textbf{,} 99-106.

Nurdiani, R. \& Zeng, C. 2007. Effects of temperature and salinity on
the survival and development of mud crab, Scylla serrata (Forsskal),
larvae. \emph{Aquaculture Research,} 38\textbf{,} 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. \emph{Fisheries Research,} 80\textbf{,} 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. \emph{Journal of Crustacean Biology,}
32\textbf{,} 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. \emph{Geochemistry, Geophysics, Geosystems,}
4\textbf{,} 1-20.

Palumbi, S. R. 2003. Population Genetics, Demographic Connectivity, and
the Design of Marine Reserves. \emph{Ecological Applications,}
13\textbf{,} 146-158.

Paradis, E. 2010. pegas: an R package for population genetics with an
integrated-modular approach. \emph{Bioinformatics,} 26\textbf{,} 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.
\emph{Environmental Modelling \& Software,} 42\textbf{,} 47-54.

Pata, P. R. \& Yniguez, A. T. 2019. Larval connectivity patterns of the
North Indo-West Pacific coral reefs. \emph{PLoS One,} 14\textbf{,}
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. \& Melia, 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. \emph{Ecology and Evolution,} 7\textbf{,} 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. \emph{PloS
one,} 7\textbf{,} e37135.

Quinn, N. \& Kojis, B. 1987. Reproductive Biology of Scylla Spp.
(Crustacea: Portunidae) from the Labu Estuary in Papua New Guinea.
\emph{Bulletin of Marine Science,} 41\textbf{,} 234-241.

Ravago-Gotanco, R. \& Kim, K. M. 2019. Regional genetic structure of
sandfish Holothuria (Metriatyla) scabra populations across the
Philippine archipelago. \emph{Fisheries Research,} 209\textbf{,}
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
\textless{}I\textgreater{}(Dascyllus aruanus)\textless{}/I\textgreater{}
in the Coral Triangle. \emph{Bulletin of Marine Science,} 90\textbf{,}
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. \emph{bioRxiv}\textbf{,}
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. \emph{Journal of Climate,} 20\textbf{,} 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. \emph{Diversity and
Distributions,} 25\textbf{,} 1684-1696.

Riginos, C. \& Liggins, L. 2013. Seascape Genetics: Populations,
Individuals, and Genes Marooned and Adrift. \emph{Geography Compass,}
7\textbf{,} 197-216.

Robertson, W. D. \& Kruger, A. 1994. Size at Maturity, Mating and
Spawning in the Portunid Crab Scylla serrata (Forskal) in Natal, South
Africa. \emph{Estuarine, Coastal and Shelf Science,} 39\textbf{,}
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. \emph{Mol Ecol,} 28\textbf{,} 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. \emph{Mol Biol Rep,} 40\textbf{,}
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.
\emph{Mol Ecol,} 27\textbf{,} 1603-1620.

Sanford, E. \& Kelly, M. W. 2011. Local adaptation in marine
invertebrates. \emph{Annual review of marine science,} 3\textbf{,}
509-535.

Schiavina, M., Marino, I. a. M., Zane, L. \& Melia, P. 2014. Matching
oceanography and genetics at the basin scale. Seascape connectivity of
the Mediterranean shore crab in the Adriatic Sea. \emph{Molecular
Ecology,} 23\textbf{,} 5496-5507.

Schunter, C., Carreras-Carbonell, J., Macpherson, E., Tintore, J.,
Vidal-Vijande, E., Pascual, A., Guidetti, P. \& Pascual, M. 2011.
Matching genetics with oceanography: Directional gene flow in a
Mediterranean fish species. \emph{Molecular Ecology,} 20\textbf{,}
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. \emph{Marine Ecology Progress Series,} 554\textbf{,} 1-19.

Sjoqvist, 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. \emph{Molecular Ecology,} 24\textbf{,} 2871-2885.

Swearer, S. E., Treml, E. A. \& Shima, J. S. 2019. A Review of
Biophysical Models of Marine Larval Dispersal. \emph{In:} HAWKINS, S.
J., ALLCOCK, A. L., BATES, A. E., FIRTH, L. B., SMITH, I. P., SWEARER,
S. E. \& TODD, P. A. (eds.) \emph{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. \emph{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). \emph{International Journal of
Fisheries and Aquatic Studies,} 2\textbf{,} 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. \emph{Mol Med,} 15\textbf{,}
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.
\emph{PeerJ,} 1\textbf{,} 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. \emph{Coral Reefs,} 36\textbf{,} 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.
\emph{Evol Appl,} 10\textbf{,} 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. \emph{Ecol Evol,} 8\textbf{,}
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. \emph{J Evol Biol,} 25\textbf{,}
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.
\emph{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).
\emph{PLoS One,} 8\textbf{,} e83493.

Wang, S., Meyer, E., Mckay, J. K. \& Matz, M. V. 2012. 2b-RAD: a simple
and flexible method for genome-wide genotyping. \emph{Nat Methods,}
9\textbf{,} 808-10.

Waples, R. S. 2010. Spatial-temporal stratifications in natural
populations and how they affect understanding and estimation of
effective population size. \emph{Mol Ecol Resour,} 10\textbf{,} 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. \emph{Evolution,} 38\textbf{,}
1358-1370.

Whitlock, M. C. \& Mccauley, D. E. 1999. Indirect measures of gene flow
and migration: FST not equal to 1/(4Nm + 1). \emph{Heredity (Edinb),} 82
( Pt 2)\textbf{,} 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). \emph{Mol Ecol,} 27\textbf{,} 2347-2364.

Zhang, Z., Schwartz, S., Wagner, L. \& Miller, W. 2000. A greedy
algorithm for aligning DNA sequences. \emph{J Comput Biol,} 7\textbf{,}
203-14.

~

\section*{Tables (separate file)}\label{tables-separate-file}

Table 1. Summary information of \emph{S. olivacea} sample location and
genetic diversity estimates.

~

Table 2. Filtering processes and specific parameters used to filter the
final SNP panel of \emph{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 \textless{} 0.05) are included in the final model.
\emph{R}\textsuperscript{2}\textsubscript{adj} represents the adjusted
coefficient of determination with p-values calculated using the analysis
of variance (\emph{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.

~

\section*{Figure legends}\label{figure-legends}

Figure 1. Estimated larval dispersal and connectivity of \emph{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 \emph{S. olivacea} connectivity based on
(a) empirical genetic estimates (pairwise \emph{F}\textsubscript{ST})
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
(\emph{F}\textsubscript{ST}) of \emph{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 \emph{F}\textsubscript{ST} 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 (\emph{p} \textless{} 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 \emph{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.

~

\section*{Data Accessibility
Statement}\label{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:
\url{https://doi.org/10.5061/dryad.3xsj3txdz}

~

\section*{Author Contribution Section}\label{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)

~

\section*{Acknowledgements}\label{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.

~

\section*{Ethical Approval}\label{ethical-approval}

\emph{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.

~

\section*{Conflict of Interest}\label{conflict-of-interest}

None declared.

\selectlanguage{english}
\FloatBarrier
\end{document}

