Adult Holothuria scabra (n=644) were sampled from natural populations at sixteen locations in the Philippine archipelago, spanning all six marine biogeographic regions (Aliño and Gomez, 1994), three island groups and fifteen Provinces ( Figure 1A ; Supplementary Table 1 ). Sites were sampled within the biogeographic regions of (1) the South China Sea: Currimao, Ilocos Norte (CUR); Anda, Pangasinan (AND); Masinloc, Zambales (MAS); Morong, Bataan (MOR); (2) Sulu Sea: Mindoro: Magsaysay (MIN); (3) Internal Seas: San Lorenzo, Guimaras (SLO); Sagay, Negros Occidental (SAG); Macrohon, Southern Leyte Gulf (SLE); Plaridel, Misamis Occidental (PLA); Tubajon, Laguindingan (TBJ); Camiguin, Misamis Oriental (CAM); (4) Celebes Sea: Olutanga, Zamboanga Sibugay (OZS); (5) North Philippine Sea: Santa Ana (STA); Laoang, Northern Samar (NSM); (6) South Philippine Sea: Guiuan, Eastern Samar (GUI) and Bislig, Surigao del Sur (SDS).

Permission to collect wild individuals, sample their tissues and transport their derivatives (DNA) was obtained by coordination with respective local governments at sites through prior informed consent, memoranda of agreements, letters of intent, or through collection and transport permits from the Philippines’ Department of Agriculture, Bureau of Fisheries and Aquatic Resources (DA-BFAR: ECC No. 2020–0094). Specimens were collected by hand from intertidal habitats and handled as outlined by Altamirano and Rodriguez (2022), being maintained in plastic storage tubs with seawater while awaiting tissue sampling. Tissue samples measuring approximately 2.0 cm × 1.0 cm were collected non-destructively, by removing a strip of skin and body wall tissue from the ventrolateral anterior flanks of specimens using a sanitized scalpel blade. Specimens were maintained out of water for no more than 2–3 minutes and returned to their collection locations following tissue sample collection, excluding individuals collected by fishers for trade. Tissues were preserved in 96% ethanol or salt-saturated DMSO buffer and maintained under refrigeration (4°C) until DNA extraction at the Marine Science Institute, University of the Philippines at Diliman, Manila, the Philippines. Total genomic DNA was extracted using a Qiagen DNEasy Blood and Tissue kit (Qiagen, Valencia, California, USA) as per Ordoñez and Ravago-Gotanco (2024), quality-checked and quantified using agarose gel electrophoresis and Qubit 2.0 Fluorometer (ThermoFisher Scientific, Wilmington, USA), respectively, before being submitted for DArTseq™ genotyping to Diversity Arrays Technology Ltd (DArT PL), in Canberra, Australia.

Diversity Arrays Technology (DArT PL) proprietary genotyping by sequencing (DArTseq™) reduced-representation libraries were prepared as described by Kilian et al. (2012) and Sansaloni et al. (2011), with a number of modifications optimised for the sandfish genome as described by Lal et al. (2021). Briefly, genome complexity reduction was achieved with a double restriction enzyme (RE) digestion, using a PstI and SphI methylation-sensitive RE combination. Custom proprietary barcoded adapters (6–9 bp) were ligated to RE cut-site overhangs as per Kilian et al. (2012), with the adapters designed to modify RE cut sites following ligation. Target “mixed” fragments (Ren et al., 2015) were selectively amplified using custom designed primers for each sample, under the following PCR conditions: initial denaturation at 94°C for 1 min, then 30 cycles of 94°C for 20 s, 58°C for 30 s and 72°C for 45 s, followed by a final extension step at 72°C for 7 min. Amplified samples were subsequently cleaned using a GenElute PCR Clean-up Kit (Sigma-Aldrich, cat.# NA1020–1KT). All samples were each normalised and pooled using an automated liquid handler (TECAN, Freedom EVO150), at equimolar ratios for sequencing on the Illumina HiSeq 2500 platform. After cluster generation and amplification (HiSeq SR Cluster Kit V4 cBOT, cat.# GD-401–4001), 77 bp single-end sequencing was performed at the DArT PL facility in Canberra, Australia.

Illumina CASAVA v.1.8.2 software was used for initial assessment of read quality, sequence representation and generation of FASTQ files from raw reads. Filtered FASTQ files were subsequently used for further filtering, variant calling and calling of final genotypes using the DArT PL proprietary software pipeline DArTtoolbox (Cruz et al., 2013). Within DArTtoolbox, the package DArTsoft14 was first used to prune reads with a quality score <25 from further processing. Individual samples were then de-multiplexed by barcode, and subsequently aligned and matched to catalogued sequences in both NCBI GenBank and DArTdb custom databases to check for viral and bacterial contamination. Any matches for viral or bacterial sequences were subsequently removed. SNP and reference allele loci were identified in reduced-representation loci (RRL) clusters and assigned scores of: “0”=reference allele homozygote, “1”=SNP allele homozygote and “2”=heterozygote. To ensure accurate variant calling, all monomorphic clusters were removed and SNP loci had to be present in both allelic states (homozygous and heterozygous). A genetic similarity matrix was produced using the first 10,000 SNPs called to assess technical replication error (Robasky et al., 2014), and furthermore clusters containing tri-allelic SNPs, aberrant SNPs and overrepresented sequences were pruned from the dataset.

Following SNP calling, each SNP was assessed for homozygote and heterozygote call rate, frequency, polymorphic information content (PIC), average SNP count, read depth and repeatability using DArT PL’s KD Compute package. These metrics were included in the raw genotype dataset supplied by DArT PL. The raw genotype dataset was further filtered to retain only a single, highly informative SNP at each genomic locus for population genomic analyses. A nested criteria of call rate and average polymorphic information content (PIC, highest to lowest rankings for both criteria) removed duplicated loci with identical Clone IDs. The dataset was further filtered for global call rate (95%), read depth (>8), average PIC (1%) MAF (2% per population) and average repeatability (95%). All loci were assessed for departure from Hardy-Weinberg Equilibrium (HWE) using Arlequin v.3.5.1.3 (Excoffier et al., 2005), using an exact test with 10,000 steps in the Markov Chain and 100,000 dememorisations. Identification of loci under potential selection was achieved using two software packages, BayeScan v.2.1 (Foll and Gaggiotti, 2008; Foll, 2012) and HacDivSel (Carvajal-Rodríguez, 2017), independently.

Bayescan computations were carried out at 0.01, 0.05, 0.1, 0.2 and 0.5 False Discovery Rate (FDR) thresholds as per Lal et al. (2021). The extreme-outlier (EOS) test in HacDivSel was used, where an algorithm searches for outlier clusters and performs a conditional LK test where more than one cluster is found based on the magnitude of F _st values. Putative F _ST outlier loci jointly-identified by both methods were removed and segregated into a separate outlier marker dataset, and all remaining SNPs tested for departure from Hardy-Weinberg Equilibrium (HWE) in Arlequin v.3.5.1.3 (Excoffier et al., 2005), using an exact test with 10,000 steps in the Markov Chain and 100,000 dememorisations. These steps collectively produced a final dataset containing selectively neutral loci.

To examine genetic relationships among populations, pairwise F _st estimates for each sample group were calculated using Arlequin v.3.5.1.3 with 10,000 permutations (Excoffier et al., 2005) and Nei’s standard genetic distances (D _S, 1978) in Genetix v.4.05.2 (Belkhir et al., 1996). A heatmap of pairwise F _st values with an associated dendrogram showing hierarchical clustering of samples was also generated using the heatmap.2 function in the R package gplot v3.1.3.1 (Warnes et al., 2020). For assessment of an isolation-by-distance (IBD) pattern of genetic structure (Wright, 1943), Mantel Tests were performed in the R package vegan (Oksanen et al., 2022), using the ‘mantel’ function which calculates a Mantel statistic based on Pearson’s product-moment correlation with 999 permutations. A parallel evaluation was carried out in the R package adegenet (Jombart and Ahmed, 2011; Jombart and Collins, 2015), using the ‘mantel.randtest’ function to compare values. Population pairwise genetic distances were computed in the R packages hierfstat (Goudet, 2005) and mmod (Winter, 2012), for F _st (Nei, 1987) and Jost’s D (Jost, 2008) measures, respectively. Geographic distances between tissue sampling sites were estimated using the path tool in Google Earth Pro v. 7.3.6.9345 (Keyhole Incorporated, 2022) and outputs of the DisperGPU particle model for dispersal pathways mediated by major surface circulation patterns (refer to section on hydrodynamic particle dispersal modelling). Correlograms of genetic vs. geographic distances were generated using the R function plot, and fitted with a LOESS-smoothed local mean line to assess relationship linearity (Figueira, 2019; Alpaydin, 2020). For Jost’s D computations, it was necessary to log-transform genetic and geographic distances to improve linearity (Wick, 2018). Local density values were generated using the 2-dimensional kernel estimation function ‘kde2d’ in the MASS R package (Venables and Ripley, 2002), and superimposed through a colour palette on correlograms to better visualise relationships.

A Discriminant Analysis of Principal Components (DAPC) was employed to investigate genetic structure using the R package adegenet 1.4.2 (Jombart, 2008; Jombart et al., 2010; Jombart and Ahmed, 2011), at country-wide (all 16 sites) and regional levels (eastern and southern vs. western coastline site subsets), utilising both selectively neutral and putatively adaptive datasets. The DAPC used an α-score optimisation to inform the number of principal components to retain and the ‘find.clusters’ function of adegenet determined the optimal number of clusters using a Bayesian Information Criterion (BIC) method (k-means clustering). As group memberships of a DAPC may be subject to overfitting the discriminant functions when too many principal components are retained, an alternative method of visualising clusters not reliant on k-means clustering was utilised (see Jombart and Collins, 2015). This permitted an additional visualisation of biologically-defined group memberships relying on retained discriminant functions only, using the ‘assignplot’ and ‘compoplot’ functions of adegenet, with the latter generating a STRUCTURE-like barplot of group membership probability.

For estimation of proportional ancestral contributions and population stratification among sample groups, the ADMIXTURE package was used (Alexander et al., 2009; Zhou et al., 2011). ADMIXTURE employs the likelihood model utilised in the STRUCTURE analysis package (Pritchard et al., 2000); however, instead of adopting a Bayesian approach and a Markov chain Monte Carlo (MCMC) algorithm to sample posterior likelihood distributions, it applies a Maximum Likelihood method to estimate parameters (Alexander et al., 2009). The optimal number of groups (K) was selected based on the K value (ranging from 2 to 18) with the lowest cross-validation error estimated across 10 independent runs (–cv=10). Individual ancestry coefficients were then calculated at the optimal k and 100 independent runs. ADMIXTURE barplots were generated using the R package pophelper v2.3.1 (Francis, 2017).

To investigate fine-scale stratification at selectively neutral loci, network analyses were also carried out using the Netview R package (Neuditschko et al., 2012; Steinig et al., 2016). Netview R population networks were generated based on a shared allele 1-identity-by-state (IBS) distance matrix created in the PLINK v.1.07 toolset (Purcell et al., 2007). Each network was constructed following computation of the maximum number of nearest neighbours for each individual (Tsafrir et al., 2005; Neuditschko et al., 2012; Steinig et al., 2016). Individual networks were then visualised and edited in the Cytoscape v.2.8.3 network construction package (Smoot et al., 2011). The IBS matrices and corresponding networks were constructed at various thresholds of the maximum number of nearest neighbour (mk-NN) values between 1 and 50, after which the optimal network for representation was selected based on cluster stability (Steinig et al., 2016).

To evaluate partitioning of genetic structure among and within populations, we carried out an analysis of molecular variance (AMOVA) in Arlequin v.3.5.1.3 (Excoffier et al., 2005) using three groupings, after the approach described by Ravago-Gotanco and Kim (2019). Using 10,000 permutations, hierarchical fixation indices and proportions of variance were examined at (1) all sample groups assessed together; (2) groups delineated by marine biogeographic regions according to respective seas and (3) groups determined by Netview R, DAPC and ADMIXTURE analyses of fine-scale genetic structure. Assessment of the direction and magnitude of migration between populations (at selectively neutral SNPs) was carried out using the ‘divMigrate’ function and Nei’s G_st method originally implemented in the now deprecated R package diveRsity (Keenan et al., 2013; Sundqvist et al., 2016). A shiny (Chang et al., 2023) implementation of divMigrate is available online as divMigrate-online at https://popgen.shinyapps.io/divMigrate-online/, and was used for analyses.

Allelic diversity indices for each sample group including the average observed (H _o) and average expected heterozygosity corrected for population sample size (H _n.b.), were computed using GenAlEx v6.502 (Peakall and Smouse, 2006) and validated with Genetix v.4.05.2 (Belkhir et al., 1996). GenAlEx was also used to determine Wright’s inbreeding coefficients (F _is), mean numbers of alleles per locus (A, MAF ≥ 5%), the effective number of alleles (N_eff) per population, number of locally common alleles (MAF ≥5%) found in <25% of all populations tested and percentages of polymorphic loci. The number of private alleles per population (A _p, at MAF ≥5%) was computed with HP-RARE v.1.0 using the rarefaction method (Kalinowski, 2004), and verified with the R package PopGenKit (Paquette, 2012).

Mean Internal relatedness (IR) within each population was calculated using the R package Rhh (Alho et al., 2010), as described by Amos et al. (2001). This metric is a modification of the Queller and Goodnight (1989) measure which examines the proportions of shared alleles, with higher weighting assigned to shared rare alleles compared to shared common alleles. Mean IR values centred around zero are suggestive of individuals descended from ‘unrelated’ parents, while negative and highly positive values suggest the presence of relatively ‘outbred’ and possibly ‘inbred’ individuals, respectively (Amos et al., 2001). Effective population sizes for each sample group based on the linkage disequilibrium method (N_eLD ) were estimated using NeEstimator v.2.01 (Do et al., 2014).

The raw DArTseq report contained a total of 85,271 SNPs genotyped across 644 individuals, and with the removal of 26,707 duplicate (clone) SNPs at each genomic locus, 58,564 SNPs were retained (31.3% reduction). Read depths ranged from 2.5–133.0 (median read depth=18.3); and further filtering for global call rate (95%), read depth (>7), average PIC (1%) MAF (2% per population) and average repeatability (95%), resulted in 10,642 SNPs remaining in the dataset. Screening for jointly-identified F _st outliers detected 117 SNPs between Bayescan (n=158 at an FDR=0.01) and HacDivSel (extreme outliers, n=122), which were separated into a dataset of putatively-adaptive loci to be used for investigation of selective signatures in genetic structure ( Supplementary Data 3 ). Testing for HWE detected a total of 103 loci deviating from expectation under HWE at p<0.00001, however none were commonly out of HWE across all populations and therefore no loci were removed from the dataset. These steps collectively produced a dataset containing 8,266 filtered selectively-neutral SNPs ( Supplementary Data 2 ), which were used for performing population genomic analyses.

The overall pattern of population differentiation at the country level for H. scabra in the Philippines is complex, with genetic connectivity highest within proximate marine biogeographic regions, and greater separation evident between geographically distant sites. Examination of broad-scale selectively neutral genetic structure among all sites with Netview R revealed separation of sites into four clusters overall ( Figure 1B ; Supplementary Figure 1 ). Individuals sampled at OZS segregated into one discrete, compact and relatively isolated cluster, while all remaining sites resolved three more diffuse groups. Two of these clusters included sandfish sampled from the North (STA, NSM) and South Philippine Seas (GUI and SDS), respectively, indicating the presence of restricted gene flow along the Eastern seaboard of the Philippines. The remaining larger cluster comprised all remaining sites from the South China, Sulu and Internal Seas (Sibuyan, Visayan and Bohol Seas) regions; CUR, AND, MAS, MOR, MIN, SLO, SAG, PLA, CAM, TBJ and SLE, suggesting high levels of gene flow and population connectivity within and among these regions. Within this cluster, sandfish sampled at CUR and AND displayed a small degree of segregation, possibly indicative of weak geographical isolation at these peninsular locations.

Examination of sites separately by marine biogeographic region provided further insights into genetic structure, with two additional pairs of DAPCs including only (1) South China, Sulu and Internal Seas sites (with the addition of STA as an outgroup), and (2) North and South Philippine, and Celebes Seas (with the addition of SLE as an outgroup). At selectively-neutral markers along the Northern and Western Philippine seaboard together with the Internal Seas region, divergence of STA in a discrete cluster and segregation of CUR, AND with MAS in a separate cluster was evident ( Figure 2D ). Concordant patterns of differentiation were apparent when individual group membership probability was examined using a DAPC compoplot ( Figure 2H ), however, at putatively adaptive SNPs, only STA remained isolated ( Figure 2E ).

Collectively at 8,266 selectively neutral markers, a total of seven genetic stocks are apparent within the three broader groups: (1) CUR, AND, MAS; (2) MOR, MIN, SLO; (3) SAG, SLE, PLA, CAM, TBJ; (4) STA, NSM; (5) GUI; (6) SDS and (7) OZS ( Figures 2B, D, F, H, I ). Furthermore, evaluation of proportional ancestral contributions through ADMIXTURE analyses resolved the optimal number of clusters as k=7 ( Figure 3 ).

Population differentiation estimates were shallow (F _st range =-0.005–0.045, mean=0.016) but significant for all but 14 out of 120 pairwise comparisons ( Table 1 ; Supplementary Figure 1 ), and were concordant with Netview R clustering patterns. Pairwise F _st estimates were highest between OZS and STA, NSM and GUI, respectively (F _st=0.043, 0.041 and 0.045, respectively); along with GUI and AND (F _st=0.040). Moderate levels of differentiation were evident between similarly geographically isolated locations including AND vs. STA and NSM, respectively (F _st=0.039 for both comparisons), SDS and OZS vs. AND, respectively (F _st=0.033 and 0.036, respectively) and OZS vs. SDS (F _st=0.031). Pairwise F _st estimates between the majority of remaining comparisons (n=68) were shallow (mean F _st=0.008, range=-0.005–0.019), suggesting high gene flow across these sites. Similar patterns of differentiation shown by pairwise F _st data were evident in Nei’s standard genetic distance (D _S) estimates between sites ( Table 1 ), with greater genetic distances between geographically distant locations. The highest estimates were again obtained between AND and STA, NSM and GUI, respectively (D _S=0.021, 0.020 and 0.020, respectively); followed by OZS against STA, AND and GUI (D _S=0.019 for all three pairwise comparisons); together with AND vs. SDS (D _S=0.018) and OZS vs. MOR (D _S=0.017).

Partitioning of genetic variance by AMOVA showed that when all sampling sites were considered together, variability among sites (F_ct ) was highest and represented 1.63% of the total genetic variance, with the remainder being distributed among individuals (98.37%, p=0.016 Table 2 ). Grouping sites into clusters determined by genetic structure analyses revealed a similar degree of partitioning, with 1.53% distributed among clusters (F_ct ), 0.49% among populations within clusters (F_sc ) and 97.99% within populations (F_st ); p ≤ 0.02. The lowest proportion of partitioning into groups was evident when sites were distributed into marine biogeographic regions segregated by sea, with just 1.27% of the genetic variance attributed (F_ct ), 0.49% among populations within biogeographic regions (F_sc ) and 98.08% within populations (F_st ); p ≤ 0.02.

Population genetic structure at putatively-adaptive markers demonstrated varying degrees of separation at regional levels ( Figure 2 ). Visualisation of genetic structure including all 16 sites in the analysis using a DAPC and 8,266 selectively-neutral SNPs reveals a similar pattern as shown by Netview R ( Figure 2B ), with strong separation of OZS into a discrete cluster, and STA, NSM, GUI and SDS in a second cluster. All remaining sites comprised a larger and compact third cluster. Visualisation using 117 F_st outlier SNPs and all 16 sites however, reveals a single, diffuse, cluster; with weaker separation of OZS individuals only ( Figure 2C ).

Along Eastern and Southern seaboard sites at neutral SNPs, OZS expectedly remained strongly differentiated, however sandfish sampled at GUI, SLE and SDS separated into their individual clusters, with NSM and STA grouping together, indicating the presence of gene flow between these two locations ( Figure 2F ). Evaluation of individual group membership probabilities supported this pattern of genetic differentiation ( Figure 2I ). Contrastingly, visualisation of genetic structure at putatively-adaptive markers across these six sites revealed separation of OZS and GUI only ( Figure 2G ). BLAST searches for functional annotations using flanking sequences associated with the F _st outlier SNPs did not return any informative results.

A moderate pattern of isolation by distance was detected by both genetic distance measures used with Mantel Tests considering all 16 populations; Nei’s (1987) F _st ( Figure 4A ) and Jost’s (2008) D ( Figure 4B ). Pearson product-moment correlation based on Jost’s D was highest when log transformed (R² = 0.382, p=0.003, Figure 4B ), compared to untransformed distances, R² = 0.307, p=0.005. Correlations were similar and significant using F _st values (R^{2 =} 0.367, p=0.002, Figure 4A ).

Gene flow patterns between and among sample sites were estimated using relative directional migration rates (m) in the divMigrate function ( Figure 5 ), and examined different hierarchical site groupings obtained from results of genetic structure analyses. Gene flow was greater among four sites along the Eastern seaboard of the Philippines (mean m=0.968 between STA, NSM, SDS and SLE; Figure 5B ), whereas analyses were not able to detect connectivity between GUI and OZS within this same grouping. Substantial gene flow was also detected among seven sites in the Internal and Sulu Seas regions (CAM, PLA, TBJ, SLE, SAG, SLO and MIN; mean m=0.780; Figure 5D ); however inclusion of four sites from the South China Sea and one site from the North Philippine Sea (STA) in the analysis reduced the degree of detectable gene flow among these twelve locations (mean m=0.562; STA, AND, CUR, MAS, MOR, CAM, MIN, SAG, SLO, TBJ, SLE, PLA; Figure 5A ). It is possible the peninsular locations of AND, CUR and MOR function as partial barriers to gene flow between these and other Western seaboard locations. Examination of AND, CUR, MOR and MAS in isolation reinforce this possibility (mean m=0.384, Figure 5C ).

Effective population sizes (N_eLD ) of sandfish varied by marine biogeographic region ( Table 3 ), with the largest variability (MOR: 116.8 [115.0–118.6] to AND: 9,341.7 [3,441.8-∞]) observed in the South China Sea region. The smallest N_e among all sites sampled was exhibited by individuals collected in the South Philippine Sea at GUI (N_eLD =23.4 [23.3–23.4]), whereas the remaining site sampled in this region (SDS) recorded an N_eLD =2,475.9 [2,254.8–2,744.8]. Continuing along the eastern Philippine seaboard, in the North Philippine Sea, the sites of STA (1,412.3 [1,276.8–1,579.6]) and NSM (3,249.4 [2,871.6–3,740.2]) generated more uniform effective population sizes, despite the relative geographic isolation of the former. Two biogeographic regions were sampled in only a single location, which were MIN in the Sulu Sea (N_eLD =592.9 [563.5–625.6]) and OZS in the Celebes Sea (N_eLD =1,514.9 [1,427.3–1,613.9]). Effective population sizes at six sites within the Internal Seas region were more uniform, ranging from SLO (983.5 [942.3–1,028.5]) to PLA (3,171.8 [2,800.7–3,655.8]).

Estimates of mean numbers of alleles per locus (A) were largely uniform among all sixteen sites, ranging from A=1.758 at GUI to A=1.804 at TBJ, with the exception of AND, which recorded the lowest value of A=1.697. Similar patterns were evident in the mean effective number of alleles per population (N_eff=1.370–1.389), again with the exception of AND where N_eff=1.363. Numbers of locally common alleles (<25% of all populations) varied by region, with a range of 0.074–0.086 detected at nine populations sampled from central regions (Eastern Visayas, Western Visayas) and Mindanao, with the exception of GUI (0.066). The lowest numbers of locally common alleles were detected in populations from Luzon and Mindoro (range=0.045–0.072), with the lowest values recorded at AND (0.045), MOR (0.058) and STA (0.062). Percentages of polymorphic loci remained above 94.5% at all sites (range=94.5–98.7% across thirteen locations) with the exception of three Luzon region sites; where values of 84.4, 91.1 and 92.4% were observed at AND, MOR and STA, respectively. No private alleles were detected in any populations.

Patterns in observed and expected heterozygosity were largely uniform (range of H_o =0.211–0.221) across all except for five sites; MAS, AND, SLO, SAG, and CAM. The lowest observed cf. expected heterozygosity was recorded at MAS (H_o =0.197, H _nb=0.242), followed by the two Western Visayan sites (SLO: H_o =0.202, H _nb=0.238 and SAG: H_{o =}0.204, H _nb=0.244). Values obtained for AND (H_o =0.209, H _nb=0.234) and CAM (H_o =0.209, H _nb=0.245) were identical for H_o . All values computed for H_o were lower than expectation, suggesting that heterozygote excess was not detected at any site. Inbreeding coefficient (F _is) values at the majority of sites (n=11) ranged from 0.115–0.142, with five sites exhibiting the highest values. These sites included MAS (F _is=0.194), SLO (F _is=0.162), SAG (F _is=0.173), SLE (F _is=0.153) and CAM (F _is=0.156). Mean internal relatedness (IR) values with the exception of AND (IR=0.047) were all highly positive (IR range=0.105–0.435, n=15), suggesting higher levels of relatedness among individuals at these locations. The comparatively low IR detected at AND is suggestive of the presence of individuals descended from unrelated parents, perhaps due to the influence of sea ranching activities from hatchery-produced juvenile releases. The highest levels of IR were detected at CAM, MAS and SAG (IR=0.435, 0.330 and 0.307, respectively), followed by intermediate levels at MIN, MOR, SLE and SLO (IR=0.263, 0.260, 0.245, 0.245, respectively).

Considered together, these data highlight small but measurable reductions in allelic diversity with higher levels of internal relatedness at insular sites (e.g. CAM, SLE, SLO, SAG, MIN, MAS and MOR).

Simulations of hydrodynamic particle dispersal revealed directional patterns of surface current movement patterns, which were primarily region-specific ( Figure 6 ; Supplementary Data 1 ). Along the eastern seaboard of the Philippines, the North Equatorial Current bifurcates into northern and southern arms, generating the Kurushio Current which travels north along the eastern Luzon coastline up towards the Babuyan Islands and Taiwan; and the Mindanao Current which turns south past Samar Island and flows along the eastern and southern coastlines of Mindanao into the Celebes Sea. Simulations demonstrate the northern bifurcation transported particles from NSM towards STA, supporting connectivity detected in genetic data between these sites. An anticlockwise gyre south of Catanduanes Island and north of the northern Samar coastline advected particles dispersed from NSM westward through the San Bernadino Strait.

Dispersal from GUI occurred in a southerly direction through the Leyte Gulf and into the Surigao Strait, supporting limited genetic connectivity with SLE. A small proportion of particles dispersed from NSM travelled southward, supporting connectivity with GUI, and to a lesser extent SDS. Particles seeded at SDS were transported southwards by the Mindanao Current past the Davao Gulf into the Celebes Sea, where they circulated in a clockwise gyre. Particle dispersal from OZS was restricted, with transportation occurring locally along coastlines, rather than into open water westwards towards Sibugay Bay and eastwards towards Igat Bay. Dispersal from the northern Mindanao sites of TBJ, CAM and PLA occurred in a westward direction into the Sulu Sea, with high levels of admixture among these locations in the Bohol Sea, while particles seeded at SAG and SLO dispersed in a south-westerly direction, via the Tañon and Guimaras Straits respectively, with only the latter advecting into the Panay Gulf and ultimately the Sulu Sea.

In the Luzon region, particles seeded at STA were transported north towards the Babuyan Islands in the Luzon Strait and further into the South China Sea by the Kuroshio Current, preventing admixture with western Luzon sites. Dispersal from CUR, AND, MAS and MOR followed a similar pattern, with advection of the majority of particles westwards into the South China Sea in several turbulent eddies. Limited connectivity was evident however, particularly between MAS and MOR, but to a lesser degree between CUR and AND, perhaps given their more peninsular locations and smaller-scale circulation patterns along the northern Luzon coastal which may limit dispersal along an otherwise contiguous coastline. Particles released at MOR and MIN were transported southward by the South China Sea throughflow current into the Sulu Sea, with admixture in an anticlockwise gyre with northern Mindanao particles was evident. Collectively, particle dispersal simulations between South China, Sulu and Internal Seas sites are concordant with levels of connectivity discovered in the genetic data, along with connectivity between the eastern seaboard sites of STA, NSM, GUI and SDS.

Shallow but significant genetic structure at neutral markers was identified across the 16 sites at which H. scabra were sampled, with populations segregating into four broad groups overall. Concordant with the findings of Ravago-Gotanco and Kim (2019), sandfish in the Philippines do not constitute a single genetic stock, but are separated by biogeographic region along the coastlines of: (1) the Celebes Sea, (2) North and (3) South Philippine Seas, (4) Sulu Sea and (5) South China Sea together with Internal Seas. While these larger groups were consistently resolved across all methods employed to examine genetic structure ( Figures 1 , 2 ), fine-scale within-biogeographic region genetic structuring was also evident at both selectively-neutral and putatively adaptive loci.

At neutral markers (8,266 SNPs), while the 16 sites sampled resolved four broad groups ( Figure 1B ), a total of seven genetic stocks are apparent within these divisions ( Figures 2B, D, F , 5): (1) CUR, AND, MAS; (2) MOR, MIN, SLO; (3) SAG, SLE, PLA, CAM, TBJ; (4) STA, NSM; (5) GUI; (6) SDS and (7) OZS. At putatively adaptive loci (117 SNPs), these groups coalesce into six units ( Figures 2C, E, G ): (1) CUR, AND, MAS; (2) MOR, MIN, SLO, SAG, PLA, CAM, TBJ; (3) STA, NSM; (4) GUI; (5) SLE, SDS and (6) OZS. Adaptive divergence between populations arises due to selective forces operating on a subset of outlier loci which increases mean fitness in their local environment, compared to other environments. Adaptive variation is also influenced by the microevolutionary forces of genetic drift, migration and mutation (Benestan, 2020). The sole discrepancy between the number of genetic units detected here at selectively-neutral (n=7) vs. putatively adaptive (n=6) loci is between SLE and SDS, and may be attributed to the close proximity of these sites, and a corresponding lack of habitat heterogeneity. At locations where adaptive divergence is evident, those sites are separated by large expanses of open ocean (e.g. OZS vs. all other sites) or are mutually geographically isolated e.g. STA vs. South China Sea sites.

Benestan (2020) proposes the definition of a fishery stock as a “demographically cohesive group of individuals of one species exploited in a specific area”, and a population as a “group of individuals that have similar demographic or genetic characteristics and thus will respond uniquely and independently to fishing”. Consideration of both selectively-neutral and putatively adaptive genetic structure is vital for delineation of conservation and management units (Nosil et al., 2009; Funk et al., 2012), to ensure accurate fishery assessment and administration (Waples and Gaggiotti, 2006; Reiss et al., 2009). Taking into account these guidelines and combining data reported here with the findings of Ravago-Gotanco and Kim (2019), who sampled additional sites in the Sulu Sea (Tawi Tawi, General Santos and the Davao Gulf) and Palawan archipelago (El Nido and Coron); it is likely that at least nine genetic stocks are present in the Philippines. The two Palawan sites of El Nido and Coron which resolved independently (Ravago-Gotanco and Kim, 2019), may be considered discrete genetic units in addition to the seven delineated in the current study.

The presence of population genetic structure in H. scabra has been established by previous studies. Nowland et al. (2017) utilised 8–10 microsatellite loci to examine wild sandfish sampled from 10 sites in Papua New Guinea (PNG) and 2 sites in Northern Australia, and detected differentiation between both Australian and all PNG sites (≥1,000 km apart), as well as between the two Australian locations (~1000 km apart). Other genetic marker types including allozymes (Uthicke and Benzie, 2001; Uthicke and Purcell, 2004), the mitochondrial CO1 gene (Uthicke et al., 2010) and simple sequence repeats (Sembiring et al., 2022) have also been used to investigate sandfish genetic structure at various Indo-Pacific sites (Brown et al., 2022). Uthicke and Benzie (2001) reported the presence of discrete genetic units between northern Australia and the Solomon Islands (~1,600 km apart), while Uthicke and Purcell (2004) also detected genetic structure between populations in Bali, Indonesia; 2 sites in northern Australia (Knocker Bay and Torres Strait); 2 sites in eastern Australia (Upstart and Hervey Bays); 9 sites in New Caledonia and 1 site in the Solomon Islands (located 500 to >4,200 km apart).

Contrastingly, Riquet et al. (2022) found no genetic structure among 4 sites (≤80 km apart) along the northwestern New Caledonian coastline at 9 microsatellite loci. In Fiji, Brown et al. (2022) reported the presence of three genetic stocks at 6,896 selectively-neutral and 186 putatively-adaptive SNPs across ~600 km, and Lal et al. (2021) at 8,725 selectively neutral SNPs reported substantial differentiation between Fiji and San Esteban, the Philippines (pairwise F_st =0.325, ~7,500 km apart). This variability in genetic structure may be explained by potential vs. realised dispersal of H. scabra larvae, as its 14–21 day pelagic larval duration (PLD) permits widespread dispersal over large distances (Altamirano and Rodriguez, 2022; Brown et al., 2024). However, complex geomorphology (as in the Philippines), habitat heterogeneity, ocean current pathways and the availability of suitable habitat for settlement all influence the degree of connectivity, and consequently population divergence and fragmentation (Ravago-Gotanco and Kim, 2019; Brown et al., 2022). Sandfish also have a biphasic settlement process, recruiting from plankton to seagrass leaves first upon completing larval development, and thereafter to sandy substrates where they become deposit feeders (Mercier et al., 2000b; Sinsona and Juinio-Meñez, 2019; Brown et al., 2024). This biological attribute is likely to contribute to the generation and maintenance of genetic structure in this species, as habitats without seagrass may not support recruitment of sandfish larvae, while population expansion may not be possible in areas with poor seagrass coverage.

The overall pattern of genetic structure for sandfish in the Philippines is influenced by geographic distance, variable ocean current-mediated gene flow, source and sink location geography and habitat heterogeneity across the archipelago. A moderate isolation by distance pattern was detected by Mantel Tests ( Figure 4 ), however gene flow among natural populations is directed by the strength and magnitude of larval dispersal via persistent surface ocean currents across coastal peninsular locations with varying rugosity (e.g. small islands and promontories). Population differentiation estimates were expectedly highest between geographically distant sites, e.g. OZS and STA, NSM, GUI (F _st range=0.041–0.045), whereas the majority of remaining pairwise comparisons e.g. Internal Seas sites reflected high levels of admixture (mean F _st=0.008, Table 1 ). Examination of gene flow patterns ( Figure 5 ) and independent hydrodynamic particle dispersal analyses ( Figure 6 ) support divisions evident in determinations of genetic structure, with high levels of admixture among Internal and Sulu Seas sites (CAM, PLA, TBJ, SLE, SAG, SLO and MIN) and between NSM and STA. Restrictions in dispersal are also apparent between the western Luzon locations of CUR, AND and MAS vs. Internal and Sulu Seas sites.

Studies of genetic structure in other marine organisms occurring in the Philippines have also reported segregation of populations according to marine biogeographic region. Two seagrasses (Cymodocea rotundata and Syringodium isoetifolium) were found to resolve three regional groups: (1) South China Sea, (2) North Philippine Sea and (3) central Philippines (Arriesgado et al., 2016; Kurokochi et al., 2016); and populations of a seagrass-associated rabbitfish (Siganus fuscescens) are segregated into four genetic units in (1) the South China Sea, (2) North Philippine Sea, (3) South Philippine Sea and (4) Sulu with Internal Seas (Ravago-Gotanco and Juinio-Meñez, 2010; Ravago-Gotanco and Kim, 2019). Within-region studies in the Philippines are few, however Mendiola and Ravago-Gotanco (2021) examined Sulu Sea populations of the mud crab Scylla olivacea and discovered weak yet significant genetic structure across 8 sites. These findings differ from results of other studies on marine broadcast spawners which report higher levels of panmixia among populations, including pearl oysters (Lal et al., 2016a), spiny lobster (Al-Breiki et al., 2018) and skipjack tuna (Anderson et al., 2020).

It is noteworthy that Ravago-Gotanco and Kim (2019) detected a pattern of isolation by distance (IBD) in their study considering 13 sites, and this investigation discovered a similar, moderate IBD trend in H. scabra across 16 sites. Ravago-Gotanco and Kim (2019) report strong population differentiation along both east-west and north-south clines, and particularly between Sulu Sea samples (El Nido and Coron) compared against other sites in the Philippines. East-west divergence in genetic structure has been detected in other marine taxa in the Philippines as summarised by Ravago-Gotanco and Kim (2019), which include the damselfish Dascyllus aruanus (Raynal et al., 2014), the giant clams Tridacna crocea, T. maxima and T. squamosa (DeBoer et al., 2014a, b), as well as the seahorse Hippocampus spinosissimus (Lourie et al., 2005). North-south divergence has also been observed in the giant clams T. crocea (Ravago-Gotanco et al., 2007), T. squamosa and T. maxima (DeBoer et al., 2014a, b), and the mottled spinefoot rabbitfish Siganus fuscescens (Ravago-Gotanco and Juinio-Meñez, 2010).

This arrangement of population genetic structure has arisen as a consequence of oceanic circulation patterns, marine basin boundaries, ecological factors and large landmasses in the Philippine archipelago which function as barriers to dispersal. As an example, the island of Luzon and opposing current directions along east vs. west coastlines account for east-west divergence between the South China and North Philippine Seas. Similarly, the presence of the Internal Seas region and bifurcating North Equatorial Current are barriers to dispersal between the eastern coastlines of the major landmasses of Mindanao and Luzon. Ravago-Gotanco and Kim (2019) mention the presence of strong genetic differentiation at a small spatial scale (F _st=0.210, ~120 km distance) between El Nido and Coron in the Palawan archipelago (Sulu Sea), which is attributed to fine-scale oceanic circulation within embayments promoting local larval retention rather than export. Similar separation ( Figures 2F, G ) was observed between several sites in this study, including NSM and GUI (F _st=0.02, ~215 km apart) which lie on opposing sides of the North Equatorial Current bifurcation at ~12–14° N latitude; and GUI and SLE (F _st=0.024, ~140 km apart) with similar influences maintaining genetic divergence.

Independent hydrodynamic particle dispersal simulation supported patterns of connectivity and larval transport pathways evident in genetic structure across the Philippines. Concordant with genetic analyses, higher levels of connectivity are apparent among central populations (e.g. Internal Seas region) relative to peripheral locations (e.g. North and South Philippine Seas). While areas of suitable contiguous habitat for sandfish (i.e., seagrass beds and sandy/muddy substrates) are widespread across the Philippine archipelago, they are separated by stretches of deep water, particularly towards island archipelagos, such as in the Sulu and Celebes Seas. The particle model utilised here (DisperGPU) is relatively rudimentary in comparison to newer toolsets such as the Connectivity Modelling System (Paris et al., 2013) and InfoMap (Bohlin et al., 2014) utilised by Deauna et al. (2021), which modelled larval retention, self-recruitment, settlement success and sink diversity in two valuable sea cucumbers (H. scabra and Stichopus cf. horrens) in the western central Philippines. While DisperGPU is useful for verification and comparison of connectivity detected in genomic data, a primary limitation is that larval competency behaviour such as homing behaviour and “survival” cannot be simulated. It is recommended that future fine-scale regional studies examining site selection for sea ranches or broodstock collection locations, utilise models incorporating competency behaviour for evaluating potential larval admixture between captive and wild populations.

Estimates of genomic diversity were largely consistent among 15 sites surveyed, with the exception of AND which exhibited a slightly lower mean number of alleles per locus (A=1.697) compared to other sites (range A=1.758–1.804, Table 3 ). This latter observation may be a result of extensive sea ranching activities in and around AND, which utilise hatchery-produced juveniles (Juinio-Meñez et al., 2017) that are now known to be less genetically diverse compared to wild conspecifics (Riquet et al., 2022; Brown et al., 2024). Within biogeographic regions, sites sampled along eastern coastlines (Eastern Visayas and Mindanao) were slightly more genetically diverse (observed heterozygosity, mean numbers of alleles, locally common alleles and percentages of polymorphic loci) compared to western sites (Luzon and Mindoro). This trend is consistent with expectations under the core-periphery hypothesis (CPH), which predicts that populations present in the core region of a species’ distribution maintain higher genetic diversity relative to peripheral locations where suitable habitat is less available, selection pressures are more intense, population fragmentation is greater and population connectivity weaker (Brussard, 1984; Eckert et al., 2008). Similar observations are reported for H. scabra by Ravago-Gotanco and Kim (2019) in their earlier study in the Philippines and Brown et al. (2022) at six sites in Fiji.

Ravago-Gotanco and Kim (2019) mention that Philippine samples exhibited substantially higher allelic richness and heterozygosity relative to sandfish sampled in northern Australia and PNG by Nowland et al. (2017), assessed using a common suite of 8 microsatellite loci. PNG and northern Australia are located at the southern limits of the known species distribution of H. scabra, relative to the central region of the coral triangle in the Pacific basin, and an overall decrease in genetic diversity along this north-south cline is consistent with CPH predictions. Erosion of genetic diversity between central vs. peripheral populations has been recorded for other marine broadcast spawners, including the black-lip pearl oyster Pinctada margaritifera (Lal et al., 2017), silver-lip pearl oyster Pinctada maxima (Lind et al., 2007), three seagrasses (Thalassia hemprichii, Enhalus acoroides and Cymodocea serrulata: Nakajima et al. (2014); Arriesgado et al. (2016); Hernawan et al. (2017)), Pacific cod Gadus macrocephalus (Fisher et al., 2022) and neon damselfish Pomacentrus coelestis (Liggins et al., 2015).

Higher levels of genetic diversity among reef-dwelling species in the coral triangle region occupying core habitats such as the Philippines and Indonesia, relative to peripheral areas (e.g. southern, northern and eastern Pacific basin margins), has been attributed to core regions functioning as stronger larval sinks, while marginal areas act as sources (Ravago-Gotanco and Kim, 2019). These dynamics have been supported with biophysical dispersal model simulations of larval transport (Kool et al., 2011; Lal et al., 2017; Baetscher et al., 2019; Deauna et al., 2021), however there is a need for distribution-wide information from future investigations to better understand both historical and contemporary genomic and demographic drivers of genetic structure and diversity. This is particularly important for H. scabra, given its potential vs. realised dispersal limits and adult habitat preferences. Ravago-Gotanco and Kim (2019) mention that trends in genetic diversity observed in their study were consistent with CPH predictions, with a non-significant gradient of decreasing genetic diversity from central to peripheral populations, coupled with increased genetic differentiation and restricted gene flow. These observations are concordant with results presented here, with genetic diversity metrics highest in the central regions of the Internal Seas, and lowest in the peripheral regions of the South China Sea, Mindoro and Philippine Seas ( Table 3 ).

Effective population sizes (N_eLD ) remained >1,000 for the majority of sites sampled (with the exception of five locations), suggesting these populations are less likely to become extirpated in the immediate future as a result of inbreeding depression, as their N_e values remain above the critical threshold of N_e > 100; however an N_e ≥ 1000 is required for retention of evolutionary potential (Frankham et al., 2014; Frankham, 2015). For the five sites which recorded lower values: GUI, MOR, MIN, MAS and SLO (N_eLD =23.4, 116.8, 592.9, 543.8 and 983.5, respectively), data suggest the presence of contemporary genetic bottlenecks at these locations, in particular GUI and MOR, which will require conservation measures to arrest further N_e declines. A summary of historical collection and processing locations in the Philippines summarised by Choo (2008a) since the late 1970s outline the key locations in the South China, Sulu and Internal Seas regions. These centres of the BDM trade include sites in Pangasinan, La Union, Cagayan, Zambales, Quezon, Batangas, Cebu, Negros Occidental, Surigao del Norte, South Cotabato, Sulu, southern Palawan, Misamis Occidental, West Central Visayas, northern and southern Mindanao and the Leyte Gulf.

The development and application of genome-wide SNPs for examination of genomic diversity and structure in H. scabra has been relatively recent, however the few studies reported thus far have offered valuable insights. Assessment of a suitable next-generation sequencing approach (DArTseq) for sandfish was reported by Lal et al. (2021), and used to investigate population genomic diversity and structure in Fiji (Brown et al., 2022) and the Philippines (this study). Genome-wide SNPs were also employed to examine loss of genetic diversity in hatchery-produced sandfish (Brown et al., 2024), examine growth rate variability (Ordoñez et al., 2021) and assign populations of origin for traceability studies (Ordoñez and Ravago-Gotanco, 2024).

The incorporation of genomic data into fishery resource management plans provides vital information on stock health, but particularly the number and spatial extent of stocks, the amount of genetic diversity present in populations, connectivity between neighbouring stocks, and the presence of local adaptation (Sweijd et al., 2000; Baetscher et al., 2019; Benestan, 2020). For heavily exploited taxa (such as sea cucumbers), genetic information can guide restocking efforts to recover depleted populations, and inform mariculture efforts for hatchery production (Brown et al., 2024).