High gene flow in the silverlip pearl oyster Pinctada maxima between inshore and offshore sites near Eighty Mile Beach in Western Australia

Introduction

Sustainable fisheries management relies heavily on an accurate understanding of stock-recruitment dynamics and patterns of dispersal in order to set biologically relevant spatial boundaries and harvest limits (Leis, Herwerden & Patterson, 2011; Palumbi, 2003). Most fisheries-targeted invertebrate marine species have a bipartite life cycle which consists of a relatively sessile adult stage and a dispersing larval stage, with the larval stage the primary vector for dispersal from natal sites. Larval dispersal is governed by the complex interplay between larval biology (e.g. swimming capacity, larval duration) and oceanographic process that act to restrict or promote dispersal and facilitate connectivity (Kinlan, Gaines & Lester, 2005; Levin, 2006; Weersing & Toonen, 2009). As a result, most marine populations exist along a continuum of connectivity, from highly connected populations that regularly exchange larvae, to isolated populations cut-off from the broader metapopulation (Waples & Gaggiotti, 2006). The strength of connections among local populations have broad implications for population maintenance and replenishment following intensive harvesting or severe environmental disturbances (Bernhardt & Leslie, 2013; Harrison et al., 2020).

The silverlip pearl oyster, Pinctada maxima, is a bivalve mollusc widely distributed throughout the Indo-Pacific and is the main species targeted for cultured pearl and pearl shell production in the region. Pinctada maxima is the largest of the Pinctada species used for pearl production (Rose & Baker, 1994), and the only species used for culturing pearls in Australia. Pearl production in Australia relies heavily on wild-caught individuals, and the Australian fishery represents the last remaining wild-capture pearl oyster fishery in the world. Each year, approximately 500,000 individuals are harvested from the wild, with the majority of fishing focused along Eighty Mile Beach in Western Australia (Zone 2-Fisheries Report 2016, Fisheries Research and Development Corporation, 2016; Fig. 1). Wild oysters are collected by divers on surface-supplied air and transported to commercial aquaculture farms for pearl production. Commercial diving for pearl oysters is restricted to depths less than 35 m and occurs predominantly in nearshore habitats in 8–15 m depth.

Pinctada maxima, along with other Pinctada species, is a broadcast spawning hermaphrodite with a pelagic larval duration of approximately 3–4 weeks (Gervis & Sims, 1992). This moderate larval duration in this group of commercially important oysters translates to gene flow within-regions, but seemingly insufficient to prevent regionally structured genetic populations (Arnaud-Haond et al., 2008; Lal et al., 2017; Lemer & Planes, 2014; Reisser et al., 2020; Rozenfeld et al., 2008). Nonetheless, amidst a backdrop of gene flow within regions, local populations may be structured at fine-spatial scales owing to local environment forces and the stochastic nature of recruitment in highly-fecund invertebrates (Lemer & Planes, 2014; Reisser et al., 2020). For example, complex circulation patterns associated with closed lagoon habitats in French Polynesia act to restrict dispersal and generate fine-scale genetic patchiness in the black-lipped pearl oyster P. margaritifera that does not reflect an expected pattern of isolation by distance (Lemer & Planes, 2014).

There is only limited evidence of connectivity in P. maxima between inshore and offshore sites along Eighty Mile Beach. At the regional scale, P. maxima populations are characterized by strong genetic structure between Indonesia and Australia, and show a strong decline in genetic diversity with latitude (Benzie & Smith-Keune, 2006; Johnson & Joll, 1992; Lind et al., 2007). Within Australia, strong genetic subdivisions occur between populations in Northern Territory, Queensland, and Western Australia (Benzie & Smith-Keune, 2006; Johnson & Joll, 1992) indicating an absence of panmixia at the National scale. Western Australian P. maxima populations are considered a single stock, with the exception of the most southerly populations within the Exmouth Gulf (Benzie & Smith-Keune, 2006; Johnson & Joll, 1992). In the Eighty Mile Beach area, oceanographic models predict that large tidal currents (up to 8 m tides) move larvae back and forth across the shelf and connect inshore and offshore fishing grounds (Condie & Hart, 2006). An absence of any genetic structure in the Eighty Mile Beach region supports extensive larval exchange between shallow and deep populations (Benzie & Smith-Keune, 2006); however, previous genetic studies have been limited in both sampling design and resolution of genetic markers to unravel any fine-scale patterns in population structure that may occur in this region.

Here, we re-visit the genetic structure in P. maxima from inshore and offshore sites near Eighty Mile Beach in Western Australia. We used a genotyping by sequencing approach and an extensive nested sampling design to explore patterns of connectivity between deep and shallow populations in a region of high economic and cultural importance. Based on previous genetic and oceanographic studies that included samples from this region, we hypothesize that gene flow is high between deep and shallow sites, and that P. maxima populations are not genetically structured across the Eighty Mile Beach region.

Methods

Results

Approximately 2M reads per sample (+/− 13,496 SE) were used for variant calling in DArTsoft14, which returned 11,280 binary SNPs (Table 1). After filtering for call rate, minor allele frequency, coverage, and HWE, we were left with 2,986 loci called across 664 individual oysters collected from 31 sites from 10 locations and two depths (Table S1; Figs. S1–S3). Our F_ST-based outlier scans revealed no locus to be a significant outlier and so we used all 2,986 loci for analyses of population genetic structure (Fig. S4). Despite our filtering criteria that removed loci out of HWE (p.adj < 0.05), we observed a strong and significant (p < 0.001) heterozygote deficit in the total population, with only 18% of all loci that passed QC (559 SNPs) to be in strict HWE before adjusting for multiple comparisons at the 0.05 significance level (Fig. 2A-orange points). Downstream analyses were carried out using both datasets, which produced consistent results, unless otherwise noted. These datasets and accompanying scripts and raw fastq files are available via the Open Science Framework (https://osf.io/v9u38/).

Consistent with the strong heterozygote deficit in the population, global inbreeding coefficients were high (F_IS = 0.211) and significantly different from zero at all sample sites (Table S2). Although oysters from offshore sites were generally larger in size than inshore sites (Fig. 2C; P < 0.001), we did not detect any differences in genetic diversity (as H_E) between depths (Fig. 2D). A global estimate of genetic differentiation (Nei, 1972) was low but significantly different from zero (F_ST = 0.006, CI [0.005–0.007]), ruling out the null hypothesis of panmixia. Despite a significant global F_ST value, only 3% of pairwise comparisons among sites (n = 12) were significant after correcting for multiple comparisons (Fig. 3A). More than half of these significant comparisons involved a Cape Bossut site (n = 8), which was the most northerly of our inshore sample locations, or an Offshore four site (n = 7), which was the most southerly of our offshore sample sites. None of these comparisons, however, remained significant when restricting our analyses to the subset of loci in strict HWE. Similarly, hierarchical analysis of molecular variance (AMOVA) indicated that a small but significant portion of genetic variation was attributed to differences between sites within location (0.074%; P = 0.01) and depth (0.026%; P = 0.01, Table S2); however, depth did not account for a significant portion of the genetic variation when restricting our analyses to the subset of loci in strict HWE (Table S3).

Consistent with low levels of differentiation under the AMOVA framework, there were no obvious spatial clustering of samples by depth based on PCA (Fig. 3B), all samples were admixed across sites, locations, and depth. Population-level clustering based on Nei’s (1972) genetic distance revealed a more complex relationship among sites and showed some evidence of divergence in most offshore sites, but without strong bootstrap support across all nodes (<50%; Fig. S5). Bayesian clustering in structure revealed that under the ∆K method, the most likely number of genetic clusters in our dataset was K = 2 (Fig. S6), but admixture plots showed all samples to be completely admixed among sites, locations and depths (Fig. S7). These patterns were confirmed by the BIC method that identified K = 1 as the optimal number of clusters in the data, with no spatial clustering of sites using DAPC (Fig. S8). Despite an absence of any clear spatial clustering among sites, we identified a weak and marginally significant pattern of isolation by distance (Mantel’s R = 0.119, P = 0.045; Fig. 3C). This pattern remained significant when focussing exclusively on inshore (Mantel’s R = 0.203, P = 0.013) but not offshore (Mantel’s R = 0.165, P = 0.110) sites. Finally, spatial autocorrelation analysis showed no significant positive autocorrelation (greater than random genetic similarity) at most spatial scales in the dataset (Fig. 3D), reconfirming the lack of spatial genetic structure among our sample sites.

Discussion

The overall pattern of genetic structure in the silverlip pearl oyster Pinctada maxima from inshore and offshore sites near Eighty Mile Beach in Western Australia reflects one of high gene flow. Levels of genetic differentiation among sites were low and Bayesian clustering analyses indicated that all samples formed a single genetic population with no obvious separation of inshore and offshore sites. These patterns of genetic homogeneity across out study site are consistent with Pinctada species from other regions across similar spatial scales (Arnaud-Haond et al., 2008; Lal et al., 2016; Lemer & Planes, 2014; Takeuchi et al., 2020), and of dispersing marine invertebrates more generally, that are often characterized by strong regional genetic structure layered on top of high (but stochastic) gene flow among local populations (Weersing & Toonen, 2009).

The patterns of high gene flow identified in this study are not surprising considering the limited spatial extent of our study, the moderate larval duration of the species (Rose & Baker, 1994), and the large tidal currents of the region (Condie & Hart, 2006). The continental shelf in northwest Australia extends 100s of km, resulting in some of the largest tides in the world (Lowe et al., 2015). These tides drive strong cross-shelf currents that can carry P. maxima larvae tens of kilometres (Condie & Hart, 2006). Particle dispersal modelling and spat surveys indicate that inshore populations along Eighty Mile Beach are largely self-seeding, with spawning and recruitment concentrated along the 8–15 m contour, and with intermittent dispersal offshore (Condie & Hart, 2006). Despite these findings, there is a long-standing hypothesis in the pearl industry that offshore populations represent a brood stock of larvae for inshore areas. We were not able to identify any clear directionality to gene flow between depths; however, our data clearly show that there is sufficient cross-shelf dispersal to homogenize the genetic structure of P. maxima in the region. Oysters collected from offshore sites were larger in size than inshore oysters, but these differences may reflect our different sampling strategies rather than a result of the gauntlet fishing strategy of the industry. Pearl oyster collection inshore was conducted during normal commercial diving operations, and it is possible that divers preferentially selected smaller shells for our study that were not size for culture and retained larger animals to be sold to the boat for profit. The collections offshore were only for scientific purposes and therefore size was not confounded by fishing preferences.

Our initial screening for deviations from HWE based on an adjusted significance threshold revealed a strong heterozygote deficit across thousands of loci. Heterozygote deficits are common in marine bivalves (David et al., 1997; Mcgoldrick & Hedgecock, 1997), and have been previously reported in P. maxima in northwest Australia using allozymes and microsatellite markers (Benzie & Smith-Keune, 2006; Johnson & Joll, 1992; Lind et al., 2007). An excess of homozygotes is generally attributed to the Wahlund effect, where genetically discrete populations are sampled as one (Pusack et al., 2014), or due to allele drop out, where an allele at a locus fails to amplify. Although some pairwise comparisons were significant, the low levels of genetic differentiation among sample sites rule out the influence of the Wahlund effect as the driving force behind the heterozygosity deficits observed here. Further our smallest spatial scale of sampling (hundreds of metres) was likely well below the spatial scale of genetic populations. If we rule out the possibility of Wahlund effect or allele drop out resulting in this signal, then the heterozygote deficits observed in P. maxima near Eighty Mile Beach could be the result of inbreeding and assortative mating driven by subtle differences in the timing of spawning of individuals and/or patchiness in larval recruitment (Mcgoldrick & Hedgecock, 1997). However, this is not known to occur in Pinctada, and considering the strength of the inbreeding signal, it is possible that allele dropout is causing the observed deficits in heterozygosity observed across most of the loci.

Estimates of genetic differentiation among sample sites showed that approximately 3% of pairwise comparisons were significant, most of which included a Cape Bossut site or an Offshore 4 site. It may be the case that complex currents and oceanographic features that form around the Cape act to entrap larvae and isolate nearby populations. Pearl oyster habitat in the Cape Bossut area is likely seeded by larvae to the north of Eighty Mile Beach region, possibly explaining the patterns of genetic differentiation observed here (Condie & Hart, 2006). These patterns of differentiation, however, were not consistent across replicate sites within our Cape Bossut sampling location, with one of the more northern Cape Bossut sites showing a lack of any significant genetic differentiation with other sampling sites. Interestingly, many of the significant pair-wise tests also included sites from Offshore 4–and localised current patterns would seem unlikely to explain this pattern. Notably, all significant pairwise values, however, became non-significant when we recalculated F_ST based on a small subset of loci in strict HWE, cautioning interpretations of the significant genetic differentiation observed here. Nevertheless, populations of highly fecund marine invertebrates often reflect a ‘sweepstakes’ chance of reproductive success, where a small group of individuals can account for a large proportion of the successful recruits to a given area (Christie et al., 2010). This results in patterns of chaotic genetic patchiness, where one generation of genetic patches cannot predict the following generation (Johnson & Black, 1982, 1984). It is likely that in a dynamic environment like Eighty Mile Beach, population structure is fluid and can change from year to year, and although significant differentiation arises in the population, successive dispersal and recruitment events over generations homogenize the population. As a result, the genetic patches are not discrete populations, but rather a small part of a larger genetic mosaic that fluctuates through time. Only through long-term temporal monitoring can we begin to unravel the complex interplay between larval biology, ocean currents, and the space-time continuum.

Supplemental Information