Genome-wide SNP analyses reveal high gene flow and signatures of local adaptation among the scalloped spiny lobster (Panulirus homarus) along the Omani coastline

Background The scalloped spiny lobster (Panulirus homarus) is a popular seafood commodity worldwide and an important export item from Oman. Annual catches in commercial fisheries are in serious decline, which has resulted in calls for the development of an integrated stock management approach. In Oman, the scalloped spiny lobster is currently treated as a single management unit (MU) or stock and there is an absence of information on the genetic population structure of the species that can inform management decisions, particularly at a fine-scale level. This work is the first to identify genome-wide single nucleotide polymorphisms (SNPs) for P. homarus using Diversity Arrays Technology sequencing (DArT-seq) and to elucidate any stock structure in the species. Results After stringent filtering, 7988 high utility SNPs were discovered and used to assess the genetic diversity, connectivity and structure of P. homarus populations from Al Ashkharah, Masirah Island, Duqm, Ras Madrakah, Haitam, Ashuwaymiyah, Mirbat and Dhalkut landing sites. Pairwise FST estimates revealed low differentiation among populations (pairwise FST range = − 0.0008 - 0.0021). Analysis of genetic variation using putatively directional FST outliers (504 SNPs) revealed higher and significant pairwise differentiation (p < 0.01) for all locations, with Ashuwaymiyah being the most diverged population (Ashuwaymiyah pairwise FST range = 0.0288–0.0736). Analysis of population structure using Discriminant Analysis of Principal Components (DAPC) revealed a broad admixture among P. homarus, however, Ashuwaymiyah stock appeared to be potentially under local adaptive pressures. Fine scale analysis using Netview R provided further support for the general admixture of P. homarus. Conclusions Findings here suggested that stocks of P. homarus along the Omani coastline are admixed. Yet, fishery managers need to treat the lobster stock from Ashuwaymiyah with caution as it might be subject to local adaptive pressures. We emphasize further study with larger number of samples to confirm the genetic status of the Ashuwaymiyah stock. The approach utilised in this study has high transferability in conservation and management of other marine stocks with similar biological and ecological attributes. Electronic supplementary material The online version of this article (10.1186/s12864-018-5044-8) contains supplementary material, which is available to authorized users.


Background
Severe decline of many commercial fish stocks highlights the emerging need for sustainable management plans for regulation and conservation of marine biodiversity. Managing marine stocks sustainably is a dynamic process and requires an in-depth understanding of the stock and its spatial boundaries, along with biological, ecological, evolutionary, economic, social or even political factors that influence the fishery [1,2]. While traditional fishery management plans rely on morphological and demographic aspects of a population such as growth, size, and mortality rates [3,4], appropriate management should also consider evolutionary criteria, including conservation of genetic diversity and maintenance of sustainable spawning stock biomass [2]. Recent studies have shown the complexity in the population genetic structure of many marine species [5][6][7]. Generally, marine organisms possess high genetic diversity and show weak population differentiation due to highly dispersive larval stages and relative absence of barriers to dispersal in the marine environment [8][9][10]. However, seascape factors (e.g. water currents, seafloor features and bathymetry) and environmental attributes can significantly influence rates of gene flow, connectivity and genetic structure in some species. Further, evolutionary processes like genetic drift and selection [11][12][13] continuously shape the genomes of marine organisms. For these reasons, defining the population structure of such organisms is challenging, but important for their conservation and management [14][15][16]. Current progress in the fields of genomics and computational biology doubtlessly offers a versatile platform for fishery managers to answer questions and issues related to population structure, stock boundaries and the level of divergence of marine organisms [1,16,17]. Recent reports support the successful application of genomic approaches to identify conservation or management units (MUs) of marine species [8,18,19]. Many of these utilise advanced genomic approaches, using high-throughput genotyping technologies i.e. Next generation sequencing (e.g. Illumina HiSeq and MiSeq platforms) and third generation sequencing (e.g. PacBio and Nanopore technology), to isolate large number of genetic markers suitable for inference of population differentiation and structure [20][21][22][23]. These technologies have enabled the development of panels of SNP markers to investigate interspecific hybrids [24,25], to assign individuals to populations, or to identify MUs [8,11,26]. Harnessing of genomic wide SNPs in the assessment of commercial marine stocks is a successful approach that could address many questions related to the level of genetic diversity and stability of these stocks [27].
The scalloped spiny lobster P. homarus (Linnaeus 1758) is characterised by a relatively long pelagic larval duration (PLD) of about 4.5-6.5 months [28], during which the larvae is exposed to oceanic dispersal as a result of currents and wind-shear, before metamorphosing into the puerulus stage and continuing its life as a benthic organism [28,29]. The species is distributed throughout the Indo-Pacific [30] and in the region supports valuable fisheries of considerable socio-economic importance. There are major concerns about the future of spiny lobster fisheries owing to a general decline in catch [31,32], emphasizing the need for serious efforts towards sustainable fishery management and regulation of the species. It is essential to introduce comprehensive fishery management guidelines for the species, considering a wide range of biological aspects e.g. demographic interactions of individuals and genetic structuring. Characterisation of stock boundaries and identification of population divergence will greatly support managers in deciding whether two populations should be managed together, or as separate stocks [27,33]. Many recent works of P. homarus have studied sub-species resolution, phylogeography throughout its wide range [34][35][36][37] or dispersal capabilities [38]. However, the fine-scale genetic structure of the species in many regions is still remains unrevealed [38].
Commercial spiny lobster fisheries have a long tradition in Oman, with the country currently being one of the major suppliers to the global market [39]. Of concern, however, is the observation that the annual harvest of Omani lobsters has declined dramatically from 2000 tons/year in the 1980's, to less than 485 tons in 2016 [40]. Presently, the lobster fishery management in Oman is primarily based on data from growth, mortality and catch rates and aims to increase population densities [39,41]. Despite regulations implemented by the government, they are not regularly reinforced (i.e. high incidence of illegal catch) and no clear legislation system against illegal practices [42]. There is a lack of knowledge surrounding the population genetic structure of P. homarus in Oman, its levels of genetic fitness and relatedness in this region. Hence, the lobster stock along the Omani coastline is currently treated as one single management unit. This study is the first to assess the genetic structure of P. homarus in Oman using high resolution genome-wide SNPs genotyping. The findings provide valuable insight into the connectivity of the Omani P. homarus population and will aid in the identification of management units for the fishery of this commercially important crustacean.

SNPs quality control and filtering
A primary dataset of 48,140 SNPs was filtered to retain 7988 SNPs (Additional file 1) suitable for genomic analysis (Table 1). Significant deviations from Hardy Weinberg Equilibrium (HWE) were observed across all populations (p < 0.000004 after Bonferroni correction).
After the removal of these SNP loci deviating from HWE, significant skew in estimation of genetic diversity indices (F is and H o ) was observed. Thus, indicating the status of those SNPs as putative null alleles or genotypic artifacts. Additionally, nine individuals were excluded from the dataset due to poor genotyping coverage < 80%.

Population genetic diversity
Observed heterozygosity (H o ) across populations ranged from 0.1660 to 0.1840 and were generally lower than the expected heterozygosity values (H n.b ) (0.2260-0.2333, Table 2). Average individual multilocus heterozygosity (Av.MLH) revealed similar values and distribution to H o, ranging from 0.1683 to 0.1858 (for Ashuwaymiyah and Haitam, respectively) ( Table 2). Average standardized MLH (sMLH) values slightly varied across populations and ranged from 0.9302 to 1.029. Inbreeding coefficient (F is ) was significantly high across all populations, (0.2094-0.2861, Table 2). The estimated parameters of identity disequilibrium (g 2 ) slightly differentiated from zero (0.0003-0.0038). However, this differentiation was statistically significant (i.e. 95% C. I. does not overlap zero) in two locations (Haitam and Dhalkut, Table 2). Estimated effective population size (N eLD ) varied from 5507.5 for Haitam and 10,305.8 for Al Ashkharah to an infinite value for other populations.

Population differentiation and genetic structure
In general, pairwise genetic differentiation estimates (F ST ) using 7988 SNPs indicated very low levels of genetic differentiation, with average F ST = 0.0004 (±SD = 0.1843), with only seven out of 28 pairwise comparisons being statistically significant (Table 3). AMOVA indicated an absence of hierarchical genetic structure between populations (variation of 80.47% within individuals; 19.51% among individuals; 0.01% among populations and 0% among groups (AS, MA; group1; DU, RM, group2; HA, SH, group3; MI, DA; group4)). While visualization of population structure using DAPC with all 7988 SNPs revealed two admixed genetic clusters (Additional file 2), population Network analysis with Netview R displayed only one genetic cluster (Fig. 1a, b). Similarly, high admixture was observed through the NJ tree, indicating high genetic relatedness among individuals from different geographical locations (similar branch lengths among all individuals, Fig. 1c).

Putatively selective SNPs
The genomic scan using Bayescan v.2.1 identified only one SNP as a directional outlier, which was   Table 4). Visualization of DAPC revealed that Ashuwaymiyah samples represented a distinct genetic cluster, while all other samples were grouped into a second admixed genetic cluster (Fig. 3). The same pattern was also observed by population Netview at a range from k-NN =10 to 30 and visualized at k-NN = 15 ( Fig. 4).    Despite repetitive culling of SNPs deviating from HWE, notably lower, but still significantly high F is values (0.2094-0.2861) were observed. Heterozygote deficits have been observed in many marine invertebrates [43] and are believed to be caused by genotypic artifacts, null alleles, population stratification caused through the Wahlund effect, biological and behavioural traits and selection [44][45][46]. In other cases, genotype calling errors have been reported to cause deviations from HWE [47][48][49]. To overcome these in our dataset, SNP loci which had call rate (< 0.7) and those significantly deviated from HWE or had low minor allele frequencies (MAF) were removed [48]. In this study, estimations of identity disequilibrium (g 2 ) slightly differentiated from zero (0.0003-0.0038, Table 2) with no significance except for Haitam and Dhalkut, thus, the obtained high F is values for P. homarus samples in this study are less likely due to inbreeding. Additionally, average multi locus heterozygosity estimates (Av.MLH) showed a moderate to high level of genomic heterozygosity (0.1683-0.1858) compared to some marine invertebrates e.g. Hapalochlaena maculosa (0.0800-0.1720) [50], Pinctada margaritifera (0.0520-0.1030) [11] and P. maxima (0.3030-0.3110) [51]. This may support that the significantly high F is values are unlikely due to reduction in genomic heterozygosity. Multi locus heterozygosity may reflect genome wide heterozygosity especially when thousands of genomic markers are used [52,53].
It is also unlikely that the high F is values are caused by the Wahlund effect, as Netview analyses indicated high levels of admixture (Fig. 1a, b). Possibly, the observed heterozygote deficits are caused by null alleles or selection. In our study, high missingness rate (> 0.2) of SNPs was observed for some samples from Ashuwaymiyah and were excluded, as they did not match our quality criteria. In the remaining samples (with genotyping coverage ≥0.8), monomorphic loci were observed for Ashuwaymiyah only, but not other sites. The monomorphic loci observed in Ashuwaymiyah samples could be attributed to mutations in restriction enzyme sites during DArTseq genotyping causing null alleles in populations [54]. It is also possible, that these observations were due to small sample size effect from Ashuwaymiyah. A study in our group with microsatellite markers revealed divergence of Ashuawymiyah stock in a larger sample size of 20 individuals (Delghandi et al., under prep).
Generally, null alleles may indicate occurrence of genetic variation in the form of point mutations or structural mutations (i.e. insertions/deletions) and contribute to the organism's fitness and adaptation [55]. It has been reported that the frequency of missing data is correlated with the level of genetic divergence between populations [54,56]. Other studies, however, have reported that null alleles are frequently encountered in SNP datasets and generate biases in estimation of diversity indices [57].

Population differentiation and genetic structure
In our study, pairwise differentiation estimates (F ST ) analysis using 7988 SNP loci, showed low but statistically significant (p < 0.05) differentiation between Ashuwaymiyah and other locations, except for Al Ashkharah and Mirbat (Table 3). Earlier reports showed that P. homarus possesses a genetic structure at the large Indo-Pacific scale [36] and relatively small scales, within both the north west Indian Ocean and south west Indian Ocean [38,58].
Other studies demonstrated different patterns and levels of genetic structure in spiny lobsters e.g. P. ornatus, exhibit low levels of genetic structure across an expansive distribution of the Indo-Pacific [59], but not at smaller scales i.e. South East Asia [60]. Similarly, no genetic structure was observed in Hawaiian P. penicillatus [61], while at the scale of the Indian Ocean and across its Indo-Pacific range, significant structure was revealed [62,63]. These contrasting patterns are referred to, mainly, environmental (e.g. pattern of water circulation) [38], bioecological (e.g. larval retention) [64], behavioral (e.g. spawning migrations) [60], or geographical factors (e.g. habitats patchiness) [59]. In this study, the observed genetic divergence among Omani P. homarus could be caused by distinct environmental and geographical factors in the region. The Omani coast in the Arabian sea is known to be influenced by complex water circulation which varies seasonally with the Monsoon and results in a series of eddies along the coast of Oman [65,66]. Generally, eddies might act as a larval retention system [67], limiting larval dispersal and maintaining divergence in marine populations including spiny lobsters, e.g. Jasus edwardsii [68] and P. h. rubellus [38]. Additionally, fragmentation of marine habitats across the Omani coastline by sandy stretches and absence of corals and rock reefs [44] could have contributed to this observation. Similar observations have been reported for other marine organisms with larval life stage such as Corkwing wrasse across Norwegian coastline [69] and the Omani clownfish [70]. In fact, Ashuwaymiyah is a shallow bay characterized by large rocky reefs, and its extended sea shelf (about 50 km) is less affected by the Monsoon currents, even during its extremes from June to September (National Survey Authority, Ministry of Defense, Oman, unpublished local data). This unique geography of Ashuwaymiyah could have limited the gene flow among lobsters from Ashuwaymiyah and other sites. This could be a possible explanation for the putative genetic divergence of P. homarus samples from Ashuwaymiyah. Moreover, it is not surprising to observe a genetically heterogeneous stock among other admixed stocks over a relatively small geographical scale. Such observation of fine scale genetic differentiation within relatively high gene flow environment has been widely described in a variety of marine species with planktonic early life stage and in a phenomenon known as chaotic genetic patchiness [71,72]. Examples include crown-of-thorns starfish, Acanthaster planci [73], clam, Spisula ovalis, [74], the sea urchin, Strongylocentrotus purpuratus [75], marine goby, Coryphopterus personatus [76], marine goby, pulmonate limpets, Siphonaria sp. [71], bicolour damselfish, Stegastes partitus [77] and spiny lobster, P. interruptus [78]. Therefore, it is possible that the observed slight divergence of Ashuwaymiyah is due to chaotic genetic patchiness. A typical feature of chaotic genetic patchiness is being temporal therefore, repetitive sampling from Ashuwaymiyah would be useful to clarify the current status of the genetic heterogeneity in Ashuwaymiyah.

Detecting potential selective SNPs
Genome-wide scan for F ST outliers using Bayescan v2.1 identified only one directional outlier, hence further analysis of population structure was not possible using this approach. In contrast, applying Arlequin resulted in more candidate outlier loci being identified, allowing further investigation for adaptive variation. It is common to detect less candidate outliers with Bayescan as it is a conservative approach and may fail to detect relatively low signals of selection [79][80][81]. Therefore other studies have used only the frequentist approach to perform F ST outliers analysis [82].
Assessment of population structure with putatively directional SNPs using both population Network and DAPC, revealed Ashuwaymiyah to be potentially under local adaptive pressures (Figs. 3 and 4). The utilization of the whole dataset of SNPs could not capture the low levels of genetic structure, while analyses based on F ST outliers allowed detection of selective divergence and identification of possible discrete stock. Similarly, many other studies showed that the use of F ST outliers could detect adaptive variation in the absence of broader analyses based on neutral markers [83][84][85]. A possible explanation for the observed divergence is the heterogeneous environmental attributes i.e. massive rocky reefs and shallow waters in Ashuwaymiyah, which might be the driver of this differentiation. In addition to the morphological/biological differences among P. homarus stocks [42], genetic studies with microsatellite markers revealed a significant divergence of lobsters from the Dhofar governorate, including Ashuwamiyah from Al Sharqiyah and Al Wusta governorates (Delghandi et al., under prep).

Implications for fishery management
The population of P. homarus along Oman is currently considered as a single homogenous stock with a single management regulation. This study shows that the samples from Ashuwaymiyah are genetically distinct from other broadly admixed samples, albeit at low levels. An earlier study reported that stock of P. homarus in Ashuwaymiyah was significantly differentiated in body size and that lobsters reach maturity at significantly lower sizes when compared to two geographically close locations [86]. The same study suggested a need to investigate the current fishery management further and indicated that Ashuwaymiyah site might require separate management. Additionally, a recent biological study of P. homarus in Oman revealed that the stock in Dhofar governorate differs from other stocks in size, time of spawning and number of spawning peaks/year, suggesting a need to consider spatial management of P. homarus along the Omani coastline [42]. Our study delivers for the first time genetic support for possible differentiation of the Ashuwamiyah stock from other locations across the coastline of Oman, using genome-wide markers, and that this stock might need to be considered for regional management. We recommend the conduct of further studies with larger number of samples, coupled with environmental and ecological data, to aid integrated assessment studies and potential discovery of unique management units.

Conclusions
Utilisation of genome wide SNPs to study the genetic status of P. homarus stocks in Oman provided valuable insights into the genetic status of the stock. This genomic resource is the first of its kind in P. homarus and the SNP dataset obtained in this study has allowed deep characterization of the lobster population genetic diversity, connectivity and structure in Oman. This study has revealed general admixture and high connectivity of P. homarus across the Omani coastline. Additionally, the study highlighted the potential prevalence of local adaptive pressures in Ashuwaymiyah. These findings indicate the importance of considering spatially customized management strategies for P. homarus across the coastline of Oman. Further studies of the genetic status of Ashuwaymiyah and stocks from other locations in Dhofar with adequate sampling based on different temporal periods together with ecological and environmental data about Omani coastline is required before any conclusive decision on the stock structure can be inferred.

Sampling and genomic DNA extraction
The commercial P. homarus fishery in Oman is situated between Ras Al-Hadd and Dhalkut (a distance of approximately 1100 km) (Fig. 2). Samples were obtained from eight commercial landing-sites, covering most of the distribution range of the lobster in Oman (Fig. 2). All samples were euthanized and purchased from local fishermen in March 2015 during the legal fishing season. A single walking leg was excised from wild caught P. homarus (n = 172) and preserved immediately in 95% ethanol until DNA extraction. All samples were purchased from local fishermen in March 2015 during the legal fishing season. Genomic DNA was extracted from tissue samples using a modified cetyl trimethyl ammonium bromide (CTAB)/ Chloroform-Isoamyl method [87]. DNA extracts were further purified through a Sephadex G50 (GE, 2007) column prior to quantification with a Nanodrop 1000 Spectrophotometer (Thermo Scientific).

Library preparation and genotyping
Genomic DNA extracts were standardised to 50 ng/μl, and sent for sequencing and genotyping using DArTseq™ technology, with Diversity Arrays Technology, Canberra, Australia [88,89]. Library preparation was completed as described by Kilian et al. [89] and Sansaloni et al. [90] with all P. homarus DNA samples being digested using a combination of PtsI and HpaII restriction enzymes. Multiplexed reduced representation libraries were then sequenced on the Illumina HiSeq2500 platform for 77 cycles.
To call SNPs and genotype each individual, raw Illumina HiSeq2500 data was first de-multiplexed into individual samples, based on sample-specific barcode sequences. De-multiplexed samples were then assessed for overall sequence quality, with any fragments with an average Q-score of < 25 being removed from the dataset. Sequences were also compared to public databases for identification of contaminant sequences, and any non-target sequences (including bacterial and viral fragments) were removed. SNP calling was conducted using the DArTsoft14 algorithm within the KDCompute framework developed by Diversity Arrays Technology (http://www.kddart.org/kdcompute.html), with initial calling parameters and filtering methods as described in Morse et al. [50] and Lind et al. [91].

SNPs quality control and filtering
To eliminate potentially aberrant SNPs, stringent quality controls were applied using custom python scripts within the DArTQC pipeline (https://github.com/esteinig/dartQC) [92]. Initially, all duplicated sequences with > 95% similarity were identified using CD-HIT and collapsed into a single cluster, or removed [93]. Further, SNPs with a call rate < 70% and those where technical replicates did not return a repeatability value of > 95% were also removed. Additionally, individuals and SNPs with > 20% missing data and SNPs with a Minor Allele Frequency (MAF) < 0.02 were excluded using Plink v1.07 [94].
To investigate the effect of sequencing depth, F is and H o were calculated for each population at different reads depth (Average SNP Counts) thresholds (3, 5, 7 and 10) to discover the degree of potential bias caused by lower call depths. Accordingly, four subsets of SNPs were generated at these sequencing depths. To detect potential genotyping artifacts, SNPs were tested for significant deviation from Hardy-Weinberg equilibrium (HWE) using Arlequin v.3.5.2.2 [95]. Any SNP loci which significantly deviated from HWE were excluded following Bonferroni correction (p < 0.000004). To assess the impact of deviation from HWE, F is and H o were calculated before and after removal of significantly deviated SNPs.

Population genetic diversity
To estimate the genetic diversity within populations, standard allelic diversity indices including average observed heterozygosity (H o ), average expected heterozygosity corrected for population sample size (H n.b. ) and inbreeding coefficient (F is ) were calculated using Genetix v.4.05.2 [96]. Effective population size, using a linkage disequilibrium method (N eLD ) was computed with NeEstimator [97]. To examine individual genome wide diversity and individual inbreeding, multi-locus heterozygosity (MLH) and identity disequilibrium parameter (g 2 ) were calculated for all individuals using the R package inbreedR [52].

Population differentiation and genetic structure
To assess population differentiation and genetic structure, a number of different statistical approaches were conducted. The extent of pairwise population differentiation was evaluated using Weir and Cockerham's unbiased F-statistics [98] through Genetix v.4.05.2 [96]. To assess hierarchical levels of population structuring, an analysis of molecular variance (AMOVA) using Arlequin v.3.5.2.2 [99] was calculated between sampling locations based on grouping samples in four groups (AS, MA; group1; DU,RM, group2; HA, SH, group3; MI, DA; group4). The grouping criterion was based on habitat similarity i.e. abundance of corals and rock substrates. Obviously, Al Wusta coast is dominated by sandy habitat and sparse coral colonies [100,101]. Region surrounding DU has been exposed for the last seven years to major construction and industrial influence, most probably having an impact on fragmentation of marine populations. Hence, HA being not affected of these factors, was grouped with SH, which is just below the border line between Al Wusta and Dhofar governates (Fig. 2). In addition, the function find.clusters in the R package adegenet [102] was used to determine the optimal number of clusters with the Bayesian Information Criterion (BIC) method. To assess levels of differentiation between the obtained genetic clusters, Discriminant Analysis of Principal Components (DAPC) was used. DAPC was performed using the optimum number of principal components (PCs) calculated using the α-score function in adegenet [103].
Finally, a network analysis with no prior population assumptions was performed to assess both broad and fine scale population structure using NetView R [104]. Net-View was run through the R implementation of NetView P [104,105] at a k-NN range from 25 to 50 as determined by a k-NN selection plot. Similarly, to visualise the extent of relatedness between individuals within each population and divergence among populations, a Neighbour-Joining (NJ) tree was constructed in MEGA6 [106]. The NJ tree was constructed using 1-proportion of shared alleles (1-psa) genetic distance matrix calculated in the R package adegenet using propShared function [102].

Identifying potential selective SNPs
To detect possible signatures of directional and balancing selection, detection of putatively selective outlier SNPs was performed using an F ST approach. To minimize false positive rates in identifying SNPs under selection, two independent statistical approaches were used. A Bayesian approach was implemented in Bayescan v.2.1 [107] and a frequentist approach [108] was implemented in Arlequin v.3.5.2.2 [99]. Bayescan estimates population-specific F ST coefficients by the Bayesian method described in [109] and uses a statistical cut-off based on the mode of the posterior distribution to detect SNPs under selection [107]. Bayescan v2.1 was used with 1:10 prior odds for a neutral model and all other parameters were kept as default (20 pilot runs of 5000 iterations followed by 100,000 iterations with an additional burn-in of 50,000) [110]. Once probabilities had been calculated for each locus, they were ranked from largest to smallest. SNPs with posterior probabilities ≥0.91-1, which are categorised as strong to decisive according to the Jeffery's scale [110], were retained. In addition, the Bayescan v2.1 function, plot R.r in the R v.3.3.1 was used to control the false discovery rate (FDR) of the selective markers at FDR of 0.05. SNPs were considered as outliers if their probability was > 0.9 at FDR of 0.05.
The frequentist approach in Arlequin v.3.5.2.2 was executed under a finite island model with 200,000 simulations and 100 demes simulated [99]. SNPs were considered as outliers based on their F ST and p values. SNPs were considered as directional loci if their F ST values fell within the upper 5% quantile and p < 0.05. They were considered as balancing SNPs if their F ST values fell in the lower 5% quantile and p < 0.05.
To assess the population structure based on directional outliers, a dataset of the putative outlier SNPs was generated. Broad scale population differentiation based on this SNP dataset was examined by calculating magnitude and significance of pairwise F ST comparisons using Genetix v.4.05.2. Population structure and network were examined based on the putative outlier SNPs using DAPC and Netview respectively. DAPC was visualised after retaining the optimum number of PCs and Net-View was run at k-NN range = 10 to 30. Availability of data and materials Genotypes data are provided as Additional files.

Additional files
Authors' contributions RDB participated in the investigation design of the study, optimization of DNA extraction, analyzed and interpreted the data, drafted, revised and edited the manuscript. SRK carried out the DNA purification and optimized samples prior to library preparation, provided support in the data analysis and contributed to the writing of the manuscript. HA and MSH, extracted the DNA and provided support in the data analysis. KRZ provided distinct advices in the data quality control and analysis. DRJ participated in the investigation design, provided support in the data analysis and writing the manuscript. MAA provided support in the data analysis and reviewing the manuscript. MD developed the concept of the project and the investigation design, collected samples, provided supervisory support, contributed to interpretation of data, writing and revision of the manuscript. All authors read and approved the final manuscript.

Ethics approval
All animal samples in this study complied with the Omani law on Sea Fishing and Protection of the Marine Biological Wealth RD 53-81 and fully complied with local fisheries management and marine protected area controls. All samples were purchased from local fishermen, hence no specific permits were required for the described filed sampling as the fishermen were required to comply with local laws regarding capture. The species sampled are not endangered or protected.

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests.

Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.