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.

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-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-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-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-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-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-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.

Results

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 (Fis and Ho) 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%.

Table 1

Filtering steps and SNPs counts retained after each step

	Retained SNPs count
Initial potential SNPs	48,140
Duplicated SNPs filters	39,086
Clustered SNPs filters	32,840
Call rate ≥ 0.7	23,764
Replication average ≥ 0.95	23,549
SNPs coverage ≥80%	20,421
Reads depth ≥ 5	14,695
SNPs MAF ≥ 0.02	12,589
HWE filters	7988
Retained SNPS for genomic analysis	7988

MAF minor allele frequencies, HWE Hardy Weinberg Equilibrium

Population genetic diversity

Observed heterozygosity (Ho) across populations ranged from 0.1660 to 0.1840 and were generally lower than the expected heterozygosity values (Hn.b) (0.2260–0.2333, Table 2). Average individual multilocus heterozygosity (Av.MLH) revealed similar values and distribution to Ho, 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 (Fis) was significantly high across all populations, (0.2094–0.2861, Table 2). The estimated parameters of identity disequilibrium (g2) 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) varied from 5507.5 for Haitam and 10,305.8 for Al Ashkharah to an infinite value for other populations.

Table 2

Genetic diversity indices for P. homarus for each sampling site using 7988 SNPs

Location	n	Hn.b (±SD)	Ho (±SD)	Av. MLH (±SD)	sMLH (±SD)	g2 (95% C.I.)	Fis (p < 0.001)	NeLD (95% C.I.)
Al Ashkharah	20	0.2324 (±0.1643)	0.1789 (±0.1495)	0.1807 (±0.0027)	1.001 (±0.0147)	0.0038 (0.0000–0.0088)	0.2349	10,305 (4860 - ∞)
Masirah	29	0.2277 (±0.1635)	0.1803 (±0.1481)	0.1823 (±0.0009)	1.009 (±0.0053)	0.0004 (0.0000–0.0008)	0.2113	∞ (∞ - ∞)
Duqm	19	0.2275 (±0.1691)	0.1798 (±0.1532)	0.1817 (±0.0011)	1.005 (±0.0062)	0.0003 (−0.0001–0.0008)	0.2142	∞ (∞ - ∞)
Ras Madrakah	17	0.2260 (±0.1715)	0.1799 (±0.1576)	0.1820 (±0.0013)	1.006 (±0.0073)	0.0006 (0.0000–0.0013)	0.2094	∞ (∞ - ∞)
Haitam	20	0.2333 (±0.1630)	0.1840 (±0.1509)	0.1858 (±0.0015)	1.029 (±0.0084)	0.0010 (0.0005–0.00151)	0.2160	5507 (3457–13,521)
Ashuwaymiyah	10	0.2287 (±0.1843)	0.1660 (±0.1677)	0.1683 (±0.0025)	0.9302 (±0.0144)	0.0031 (−0.0002–0.0077)	0.2861	∞ (∞ - ∞)
Mirbat	20	0.2270 (±0.1678)	0.1780 (±0.1513)	0.1797 (±0.0024)	0.9947 (±0.0131)	0.0044 (−4.8190–0.0122)	0.2203	∞ (∞ - ∞)
Dhalkut	29	0.2273 (±0.1640)	0.1770 (±0.1437	0.1788 (±0.0017)	0.9892 (±0.0094)	0.0032 (0.0007–0.0059)	0.2248	∞ (∞ - ∞)

n number of samples, Hn.b average expected heterozygosity corrected for population sample size, Ho observed heterozygosity, Av. MLH average multi-locus heterozygosity, sMLH standard multi-locus heterozygosity, g2 identity disequilibrium parameter, Fis inbreeding coefficient, N effective population size by the linkage disequilibrium method with 95% confidence interval

Population differentiation and genetic structure

In general, pairwise genetic differentiation estimates (FST) using 7988 SNPs indicated very low levels of genetic differentiation, with average FST = 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).

Table 3

Pairwise FST values for 7988 SNPs using Genetix v.4.05.2 with permuted p-values inside brackets

	Al Ashkharah	Masirah	Duqm	Ras Madrakah	Haitam	Ashuwaymiyah	Mirbat
Masirah	−0.0006 (0.970)	–	–	–	–	–	–
Duqm	0.0001 (0.570)	0.0007 (0.198)	–	–	–	–	–
Ras Madrakah	0.0004 (0.416)	−0.0002 (0.347)	0.0007 (0.306)	–	–	–	–
Haitam	0.0003 (0.319)	0.0002 (0.477)	0.0013a (0.044)	0.0001 (0.690)	–	–	–
Ashuwaymiyah	0.0006 (0.267)	0.0002a (0.045)	0.0016a (0.041)	0.0021a (0.010)	0.0021a (0.031)	–	–
Mirbat	−0.0008 (0.931)	−0.0004 (0.901)	0.0013a (0.040)	−0.0001 (0.771)	− 0.0001 (0.720)	0.0004 (0.258)	–
Dhalkut	−0.0003 (0.485)	0.0005 (0.129)	0.0005 (0.198)	−0.0001 (0.780)	0.0007 (0.112)	0.0021a (0.043)	0.0007 (0.491)

adenotes significant comparisons

Fig. 1

Population structure of 164 individuals of P. homarus samples using 7988 SNPs. Population network was constructed using NetView R v.1.0. The network is visualised at nearest neighbor (k-NN) threshold of 25 (a) and 50 (b). Un-rooted neighbor-joining tree (c) was drawn in MEGA6 using 1-psa genetic distances, scale bar indicates genetic distance

Putatively selective SNPs

The genomic scan using Bayescan v.2.1 identified only one SNP as a directional outlier, which was characterized by Jeffery’s scale as decisive. The frequentist approach with Arlequin v.3.5.2.2 indicated out of 7988 SNPs, 504 SNPs (6.31%) potentially under divergent selection, 168 SNPs (2.10%) under balancing selection and 7316 SNPs (91.59%) likely to be neutral (Genotypes of 504 SNP loci are available as Additional file 3). Global FST for the 504 SNPs candidate outliers was 0.0423 (Table 4), more than 100 times greater than for all 7988 SNPs (0.0004). All of the pairwise estimates with these directional outliers were significant (p < 0.001), with highest differentiation between Ashuwaymiyah and Ras Madrakah (0.0736) and the lowest between Masirah and Dhalkut (0.0288) (Fig. 2; 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).

Table 4

Pairwise FST values for 504 directional SNPs using Genetix v.4.05.2

	Al Ashkharah	Masirah	Duqm	Ras Madrakah	Haitam	Ashuwaymiyah	Mirbat
Masirah	0.0366	–	–	–	–	–	–
Duqm	0.0407	0.0457	–	–	–	–	–
Ras Madrakah	0.0423	0.0678	0.0489	–	–	–	–
Haitam	0.0325	0.0406	0.0446	0.0411	–	–	–
Ashuwaymiyah	0.0662	0.0678	0.0696	0.0736	0.0664	–	–
Mirbat	0.0383	0.0322	0.0476	0.0437	0.0343	0.0706	–
Dhalkut	0.0352	0.0288	0.0397	0.0387	0.0332	0.0668	0.0309

All pairwise FST values are significant (actual permuted p-values < 0.001 for all estimates)

Fig. 2

Sampling sites of P. homarus specimens along the coastline of Oman. AS Al Ashkharah, MA Masirah Island, DU Duqm, RM Ras Madrakah, HA Haitam, SH Ashuwaymiyah, MI Mirbat and DA Dhalkut. Numbers in circles represent genotyped samples in the final analyses. The black circle represents Ras Al-Had, the northern border of the commercial fishery sites of P. homarus. The map was obtained from: (https://www.d-maps.com/carte.php?num_car=516&lang=en) and edited to highlight the sampling sites

Fig. 3

DAPC scatter plot created using eight optimum PCs out of 504 outlier SNPs across 164 P. homarus individuals in the R package adegenet. Dots represent individuals. The plot showing Ashuwaymiyah represents a distinct genetic cluster and other sampling locations represent a second admixed genetic cluster

Fig. 4

Population network of P. homarus individuals using NetView R v.1.0 based on identity by similarity (IBS) distance matrix calculated in PLINK after Steinig et al. [104]. The network was drawn at nearest neighbor (k-NN) threshold of 15, using 504 potentially directional SNPs and 164 individuals

Discussion

Assessment of genetic diversity and population structure of P. homarus along the Omani coastline is vital for its stock management and conservation. This study is the first to identify a genome-wide SNP dataset for P. homarus utilising DArTseq technology and to determine the genetic stock structure of the species in Oman. The obtained SNP dataset revealed unique insights into genetic diversity and stocks boundaries of the spiny lobster in Oman. Significant deviations from Hardy Weinberg Equilibrium (HWE) were observed across all populations, apparently caused by a heterozygote deficit as represented by significantly high Fis values (all Fis values > 0, p < 0.001). Despite repetitive culling of SNPs deviating from HWE, notably lower, but still significantly high Fis 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-46]. In other cases, genotype calling errors have been reported to cause deviations from HWE [47-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 (g2) slightly differentiated from zero (0.0003–0.0038, Table 2) with no significance except for Haitam and Dhalkut, thus, the obtained high Fis 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 Fis 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 Fis 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].

In our study, pairwise differentiation estimates (FST) 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 FST 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-81]. Therefore other studies have used only the frequentist approach to perform FST 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 FST outliers allowed detection of selective divergence and identification of possible discrete stock. Similarly, many other studies showed that the use of FST outliers could detect adaptive variation in the absence of broader analyses based on neutral markers [83-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.

Methods

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].

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, Fis and Ho 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, Fis and Ho were calculated before and after removal of significantly deviated SNPs.

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

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]. NetView 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 FST 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 FST 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 FST and p values. SNPs were considered as directional loci if their FST values fell within the upper 5% quantile and p < 0.05. They were considered as balancing SNPs if their FST 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 FST 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 NetView was run at k-NN range = 10 to 30.

Genotypic data of P. homarus for 164 individuals from Oman using 7988 SNPs. Genotypes are in genetix format. (TXT 576 kb)

A plot of the Discriminant Analysis of Principal Components (DAPC) against the discriminant function retained. The plot is indicating presence of two genetic clusters of P. homarus in Oman. The plot was generated using the most informative 13 PCs identified from all 7988 SNPs dataset across 164 P. homarus individuals in the R package adegenet. (TXT 9105 kb)

Genotypic data of P. homarus for 164 individuals from Oman using 504 putatively directional SNPs. Genotypes are in genetix format. (PDF 389 kb)