Assessment of genetic structure among Australian east coast populations of snapper Chrysophrys auratus (Sparidae)

Additional keywords: effective population size, fisheries management, microsatellite genotyping, stock structure.

Chrysophrys auratus (Forster, 1801) is a highly prized fish found in coastal waters off Australia and New Zealand. Throughout Australia the species is a key target in both commercial and recreational fisheries (Kailola et al. 1993). The Australian east coast fishery has components in New South Wales (NSW), Queensland (Qld), Victoria (Vic.) and Tasmania (Tas.) that are currently assessed and managed independently within each jurisdiction (Fowler et al. 2016). Stock assessments, based on a wide array of biological and fishery data, are regularly conducted to provide fisheries managers with the information used to inform the regulation of harvest impacts. Knowledge relating to stock structure is important for stock assessment and management because it informs the spatial units (fishery regions) over which data should be integrated and management should operate. Further, snapper along the east coast of Australia spawn in coastal waters, where their larval dispersal is influenced by the East Australian Current (EAC; Curley et al. 2013). Predicting longer-term evolutionary changes, such as those related to ocean warming or changes to currents (Pecl et al. 2014), depends on knowledge of genetic structure. If populations are in decline, knowledge of genetic connectivity can assist in identifying source or sink populations (Hauser and Carvalho 2008).

The need for clarification of genetic structure of snapper along the Australian east coast was emphasised by recent stock status assessments that indicated that, despite the considerable amount of fishery and biological information available, the stock status of east coast snapper was ‘uncertain’ (Fowler et al. 2016). This uncertainty was largely because of state-based assessments of populations in Qld suggesting they were ‘recruitment overfished’, whereas in NSW they were suggested as possibly ‘growth overfished’ but not ‘recruitment overfished’ because commercial catch rates were increasing and the sizes and ages in landings remained stable (Finn et al. 2015). The status of stocks in Tas. and eastern Vic. was also undefined because of insufficient data (Kemp et al. 2012; Fowler et al. 2016). Although the major east coast C. auratus fisheries in NSW and Qld are managed independently, confirmation of whether or not both fisheries exploit the same biological stock, and thus should be combined for fisheries assessment, is critical in order to reconcile the uncertainty in exploitation status of east coast snapper (Ferrell and Sumpton 1997; Sumpton et al. 2008; Fowler et al. 2016). Further, should the limited information available for fisheries in eastern Vic. and Tas. continue to influence the assessment of the main fisheries further north?

Rapidly evolving molecular technologies have increased in power, and reduced in cost, such that they now play an important role in defining genetic stocks of exploited fish species (Carvalho and Pitcher 2012). Stocks that are genetically different from each other are described as reproductively isolated units or regions within a species distribution (Ovenden 1990). For genetic stocks to be identified, it is necessary that interbreeding among different stocks has been negligible over very long time periods. This definition may or may not relate to ‘fisheries’ or ‘harvest’ stocks, which are described as regions within a species distribution where fishing pressure in one region does not affect the abundance of fish in other regions, implying independence of the processes of replenishment, growth and mortality, but not necessarily different genetics (Gauldie 1988). Identification of genetic stocks provides highly definitive broad spatial regions on which to focus fisheries assessment and identification of management units.

The genetic structure of east coast populations of C. auratus was first investigated in 1980 using allozyme data as part of a study examining the Australian distribution of the species (MacDonald 1980). That study failed to detect differences among the three east coast locations sampled. Sumpton et al. (2008) also used allozyme loci, but sampled the east coast more intensively at seven areas spanning a finer spatial scale between Sydney in NSW and Rockhampton in Qld (Fig. 1). The authors of that study identified a weak genetic disjunction (based largely on a difference at one locus) on the central coast of NSW, north of Sydney but south of Forster (Fig. 1). A signal of isolation by distance (IBD) was found north of the weak genetic disjunction, to the Swains Reef off Rockhampton (Sumpton et al. 2008).

Fig. 1.  Sample collection areas for the Chrysophrys auratus study, with the number of unique fish genotyped in parentheses. The weak genetic disjunction detected using allozyme data (Sumpton et al. 2008) is marked with a thick dashed line. WA, Western Australia; NT, Northern Territory; Qld, Queensland; NSW, New South Wales; Vic., Victoria; Tas., Tasmania.

Chrysophrys auratus is a long-lived species (up to 40 years; Norriss and Crisafulli 2010) that can grow up to 16 kg and 120 cm in length (MacDonald 1982; Kailola et al. 1993). The species undergoes juvenile sex inversion, with small juveniles having female reproductive tissue and a proportion of individuals becoming male near the age of maturity at (3+ years and at a reasonably large size of 25–30 cm; Francis and Pankhurst 1988; Stewart et al. 2010). These biological characteristics, combined with an intermediate fecundity (Saunders et al. 2012), mean that C. auratus is more susceptible to overfishing than faster-developing species such as sardine and anchovy, particularly if exploitation is not managed at the appropriate spatial scales relevant to population replenishment. Adult reproduction is by broadcast spawning in breeding aggregations, which, depending on the geographic area, can occur in both sheltered inshore bays and coastal waters (Kailola et al. 1993). Tag–recapture and otolith chemistry research suggest that migration of juvenile and adult snapper is highly variable among the geographic regions where studies have been conducted, ranging up to at least hundreds of kilometres (Sanders 1974; Gillanders 2002; Coutin et al. 2003; Moran et al. 2003; Parsons et al. 2003; Sumpton et al. 2003; Hamer et al. 2011; Fowler et al. 2017). Along the east coast of Australia, even if adults showed restricted movements, dispersal of the 20- to 30-day larval stage in prevailing ocean currents, particularly the dominant EAC, provides considerable scope for broad-scale north to south genetic mixing (Ferrell and Sumpton 1997). Despite having a high larval dispersal potential, C. auratus could use local coastal circulation to recruit locally, or could migrate back to its natal origins before breeding, which could create complex genetic structure. Other studies investigating the genetic structure of species with similar pelagic larval periods, namely Pristipomoides multidens (Ovenden et al. 2004), Lutjanus carponotatus (Harrison et al. 2012) and Argyrosomus japonicus (Barnes et al. 2016) have detected genetic structure suggesting low-level mixing despite potential for wide larval dispersal.

Microsatellite markers offer a high level of resolving power because of their potential for high polymorphism. They have been successfully used in New Zealand stocks of C. auratus to investigate population structure (Ashton 2013) and larval recruitment into marine reserves (Le Port et al. 2017). Ashton (2013) found that despite measuring a high level of polymorphism in the 17 microsatellite loci screened, there was little genetic differentiation to distinguish among stocks of C. auratus collected over a 900-km range in New Zealand. Although the genetic results reported by Ashton (2013) supported long distance dispersal, he also reported limited genetic mixing, over a 40-year period, between two locations situated 50 km apart. He proposed site-specific differences in mobility of New Zealand C. auratus to explain the conflicting results (Ashton 2013).

Le Port et al. (2017) used the highly polymorphic microsatellite loci to screen for relatedness among samples to assess the contribution of marine reserves to larval recruitment, finding predominantly small-scale larval dispersal distances of less than 40 km. Microsatellite markers have also been used to investigate the genetic structure of C. auratus in and around Shark Bay, Western Australia (Gardner et al. 2017). As with the New Zealand work, Gardner et al. (2017) also recorded highly polymorphic loci (only 1 of the 12 microsatellites overlapped with the New Zealand studies) and failed to detect significant genetic structure over a 350-km range, despite earlier allozyme studies suggesting that complex genetic structure may exist (MacDonald 1980). Overall, to date individual-based genetic analyses suggest that C. auratus exhibit limited larval dispersal, but lack clear genetic differences, probably due to some level of interbreeding between neighbouring populations over long time periods (an IBD structure).

Genetic markers offer an independent tool that can complement and improve traditional demographic estimators of effective population size (N_e; Luikart et al. 2010). Estimating N_e is challenging in marine species that are highly fecund and have large populations with moderate gene flow (Hare et al. 2011). Despite recent methodological developments (Waples et al. 2011), the relationship between census population size (N) and N_e is still poorly understood (Ovenden et al. 2007; Waples 2016), but is discussed in detail by Ovenden et al. (2016). Hauser et al. (2002) estimated N_e in C. auratus populations from New Zealand and found that N_e was five orders of magnitude below the adult census size (N). Several suggestions have been made regarding the factors that may contribute to very low estimates of N_e, such as high fecundity and biased reproductive success (Hauser et al. 2002; Hauser and Carvalho 2008; Coscia et al. 2016; Waples 2016). Using simulated data, Waples (2016) explored, in turn, the potential effects of increased longevity, fecundity and variance in reproductive success, as well as increased egg quality with age, and found that ‘very tiny’ N_e required some version of Hedgecock’s (1994) ‘sweepstakes’ hypothesis, whereby only a few families reproduce successfully. Waples (2016) also noted that when N_e is large (>10 000), the frequency distribution of estimates of N_e will be bimodal with either infinitely large estimates or otherwise downwardly biased values ranging in the low hundreds to low thousands. Estimating unbiased N_e for empirical data with large N requires good power in terms of sample numbers and number of loci. Although correlating N_e with N may require more biological information than is currently available for C. auratus, relative comparisons of N_e among sample areas is possible and may provide some insight into the effect of local fishing pressure on a population.

The overall objective of this study was to characterise the genetic stock structure of C. auratus along the east coast of Australia using microsatellite markers and to test the specific findings of previously investigated allozyme markers (MacDonald 1980; Sumpton et al. 2008). In addition, more comprehensive sampling than previous allozyme studies has been used to examine the findings of these previous studies, namely the IBD signature, a genetic disjunction north of Sydney and the extent of panmixia among Qld and NSW sites. The genetic stock structure of east coast populations of C. auratus determined in this study will be used to inform stock assessment models currently being developed to better inform cross-jurisdictional assessment and management of east coast C. auratus. Preliminary estimates of N_e are also made to ascertain whether existing methodology is capable of making bounded predictions for C. auratus stocks and provide a potential future proxy for changes in population size.

Fresh fin clip tissue samples of C. auratus were collected from recreational and commercial fishers who donated samples from fish that were killed during harvest between 2012 and 2016. Samples were stored in individual tubes of 100% molecular-grade ethanol. Sampling was conducted along the east Australian coastline from nine areas, spanning four states and ranging over 2000 km (Fig. 1). The sampling strategy was designed to span a wide geographic range of commercial fisheries including sites from both sides of the genetic disjunction in central NSW identified in the allozyme study of Sumpton et al. (2008; marked on Fig. 1). Most Tasmanian samples (n = 21) were sourced from waters near Devonport in the central north, but because of a paucity of samples, four additional fish from waters off Stanley in the north-west and four from waters closer to Hobart in the south-east were included. All tissue samples were transferred to the Department of Agriculture and Fisheries at the Eco-Sciences Precinct in Brisbane (Qld, Australia) for molecular analysis.

Approximately 3 mm of fin clip tissue was washed in 1 mL of MilliQ water to remove the ethanol preservative before DNA extraction. DNA was extracted using a DNeasy Tissue Kit (Qiagen, Melbourne, Vic., Australia) according to the manufacturer’s instructions into a final elution volume of 50 µL. DNA concentration was quantified using a NanoDrop 8000 spectrophotometer (Thermo Fisher Scientific, Melbourne, Vic., Australia).

This project targeted nine microsatellite loci from the larger panel that Le Port et al. (2017) used for their New Zealand C. auratus study. The 17-loci panel of A. Le Port, S. Lavery, N. Kaur, and J. C. Montgomery (unpubl. data) was reduced to exclude loci that were difficult to amplify and those Le Port et al. (2017) had identified as out of Hardy–Weinberg equilibrium (HWE). Primer sequences, annealing temperatures multiplex combinations plus the original source of each locus are detailed in Table S1, available as Supplementary material to this paper.

Rather than individually labelling each primer with a fluorescent probe, the forward primer at every locus was modified with one of four M13 sequences (Table S1). For each assay, the reaction tube contained the primer pair plus a fluorescently labelled M13 primer (with FAM, NED, VIC or PET fluorophores; Oetting et al. 1995; Schuelke 2000; Missiaggia and Grattapaglia 2006; Kirchoff et al. 2008). Thus, the nine loci were polymerase chain reaction (PCR) amplified in three multiplexed PCR reaction tubes (M1–M3) and three single locus reaction tubes (S1–S3; Table S1). A Multiplex PCR Kit (Qiagen) was used to amplify the DNA in a final volume of 6 µL. The PCR samples contained 3 µL of 2× Master Mix, 0.6 µL of 5× Q solution, varying primer concentrations (detailed in Table S1, with the labelled M13 primer concentration the same as the reverse primer) and ~20 ng of Genomic DNA template. Microsatellite PCR amplifications were performed in a Bio-Rad (Sydney, NSW, Australia) thermal cycler (DNA Engine Peltier). The DNA template and enzyme were denatured at 95°C for 15 min, followed by 37 cycles of 94°C for 30 s, 52–62°C (specific annealing temperatures are given in Table S1) for 45 s and 72°C for 90 s. To ensure consistent allele calling during genotyping, a final extension at 72°C for 45 min was used to ensure complete extension of the PCR products. Allele sizing was determined using GeneScan 500 LIZ dye size standard (Thermo Fisher Scientific). Products were separated by electrophoresis on an ABI3130xl sequencer (Applied Biosystems, Thermo Fisher Scientific). According to the manufacturer’s recommendation, LIZ peaks at 35 and 250 were excluded before fragment analysis because of their temperature sensitivity; then, microsatellite peaks were scored using Geneious (ver. 8.1.9, see http://www.geneious.com, accessed 16 May 2016; Kearse et al. 2012). A repeat positive control sample was run on every 96-sample plate to ensure scoring consistency was maintained between electrophoresis runs. Samples returning low (<200 fluorescence units) or no signal strength for a subset of loci were initially subjected to another PCR run with increased starting DNA. If the signal continued to be weak, the multiplex was broken into single-locus reactions.

To summarise genetic diversity, several different metrics were determined for the loci. A relatedness screen (in GenAlEx, ver. 6.5, see http://biology-assets.anu.edu.au/GenAlEx/Welcome.html, accessed 17 November 2016; Peakall and Smouse 2006) was used to identify duplicate samples, which were removed from subsequent screening. To estimate the level of genetic differentiation between sampling locations the standardised measure of genetic differentiation, Jost’s estimate of differentiation (D_EST; Jost 2008) for each location pair was determined (using GenAlEx, ver. 6.5, see http://biology-assets.anu.edu.au/GenAlEx/Welcome.html; Peakall and Smouse 2006). D_EST provides a more accurate measure than the fixation index (F_ST), which tends to be downwardly biased by high allelic diversity at microsatellite loci, small numbers of sampling locations and low sample numbers. The number of alleles was determined to estimate the polymorphism information content (PIC) of the loci, as well as to calculate observed and expected heterozygosity values (using Cervus, ver. 3.0.7, see http://www.fieldgenetics.com/pages/aboutCervus_Overview.jsp, accessed 11 October 2016 Kalinowski et al. 2007).

To determine deviations from HWE, exact tests were used to test each locus, in each sampling location (using Genepop-on-the-Web, ver. 4.2; http://genepop.curtin.edu.au/, accessed 28 November 2017; Raymond and Rousset 1995; Rousset 2008). Although each locus by sampling location test was independent, the same null hypothesis was tested multiple times; thus, a subjective decision was made (following Cabin and Mitchell 2000) to apply an intermediate Bonferroni-type correction (set at 9) for multiple tests (Rice 1989). For loci out of HWE, the direction of bias was determined and an assessment was made to determine whether the result could be attributed to scoring errors, allele dropout or null alleles (using Microchecker, ver. 2.2.3, see http://www.nrp.ac.uk/nrp-strategic-alliances/elsa/software/microchecker/, accessed 5 October 2016; Van Oosterhout et al. 2004). Linkage disequilibrium (LD) was tested using log-likelihood ratio statistics (G-tests) to assess each pair of loci within each sampling location (using Genepop-on-the-Web, ver. 4.2, see http://genepop.curtin.edu.au/; Raymond and Rousset 1995; Rousset 2008) with Bonferroni correction. Two assumptions made by F-statistics and structure analyses are that populations are in HWE and loci are not linked. Loci failing to comply with these assumptions were removed before analysing the data.

A power analysis to assess the resolving power of the microsatellites to detect genetic differentiation was conducted in POWSIM (ver. 4.1, see http://internt.zoologi.su.se/~ryman/, accessed 5 February 2018; Ryman and Palm 2006; Gardner et al. 2015). As part of the POWSIM evaluation, N_e was set to 10 000 following Gardner et al. (2017) e this value is not based on the linkage disequilibrium based N_e (LDN_e) empirical estimates below). Time since divergence (t) was varied to obtain seven divergence levels (F_ST) between 0.0002 and 0.005. After drift, the base population was subdivided into nine populations for simulations, with the size of each population following the sample data in Fig. 1. The mean of 500 replicates was used to estimate the proportion of samples for which the F_ST values were significantly different from zero using Fisher’s exact tests.

The population genetic structure of species across the nine sampled locations in eastern Australia was investigated using four approaches. First, Bayesian inference was used to assign individuals to expected stocks using Structure (ver. 2.3.4, see https://web.stanford.edu/group/pritchardlab/structure.html, accessed 6 December 2013; Pritchard et al. 2000). The most likely number of genetic clusters was determined following Evanno et al. (2005), as described below. Second, population pairwise D_EST (GenAlEx, ver. 6.5, see http://biology-assets.anu.edu.au/GenAlEx/Welcome.html; Peakall and Smouse 2006) values were estimated for all pairs of sampling locations. Neighbouring locations with non-significant fixation values were then pooled and pairwise fixation values recalculated in an iterative approach to identify possible spatial boundaries to genetic stocks (Broderick et al. 2011). Third, a discriminant analysis of principal components (DAPC; Jombart et al. 2010) available in the Adegenet package (ver. 2.1.1, see https://cran.r-project.org/web/packages/adegenet/index.html, accessed 5 February 2018; Jombart 2008), run through RStudio (ver. 0.99.903, RStudio, Boston, MA, USA, see http://www.RStudio.com/ide, accessed 5 February 2018) was used to distinguish genetic clusters. Finally, an analysis of molecular variance (AMOVA; Arlequin, ver. 3.5.1.2, see http://cmpg.unibe.ch/software/arlequin35/, accessed 18 December 2013; Excoffier and Lischer 2010) was used to determine the percentage genetic variance explained by the groupings deduced from the structure and D_EST analyses.

For the structure analysis, a series of Markov chain Monte Carlo (MCMC) simulations was run using models of both admixture and no admixture, and using locations as priors correlated with allele frequencies (Falush et al. 2003). Simulations were run for a range of stock sizes (K = 1–10) to determine the optimal number of clusters following Evanno et al. (2005). Ten repetitions were run for each stock size, burn-in was set to 10⁴ and 10⁶ repetitions were run after burn-in. Using the ΔK estimator approach of Evanno et al. (2005), the rate of change in the log probability of the data between successive K values was calculated (ΔK) and plotted against K to determine the most likely number of genetic stocks.

Multivariate DAPC was conducted using the complete dataset because discriminate analyses are robust to loci out of HWE or in LD. Evidence of genetic clusters was examined in DAPC by running successive K-means clustering in the find.clusters function with scaling activated during the principal component analysis (PCA) to give higher influence in the clustering to loci with more alleles. The optimal number of clusters was determined as the K with the lowest Bayesian information criterion (BIC; Jombart et al. 2010). DAPC was also run for a priori stock number K = 2 based on the outcomes of the structure analysis and population pairwise D_EST analysis.

Spatial patterns of genetic divergence were investigated using a genetic model of IBD correlating genetic distance (D_EST/(1 – D_EST)) to geographic coast distance (km) in Genepop-on-the-Web (ver. 4.2, see http://genepop.curtin.edu.au/; Raymond and Rousset 1995; Rousset 2008). Shoreline distances (km) between sampling locations were estimated manually using Google Maps (API, ver. 3.30, see https://www.google.com.au/maps/ accessed 18 November 2017, Google, Mountain View, CA, USA), factoring in land barriers. The resulting correlation for all sampling locations was plotted in Microsoft Excel 2013 (ver. 15.0.5075.1000, Microsoft, Bellevue, WA, USA) and was assessed using Mantel (1967) tests and distance-based redundancy analysis (dbRDA; Legendre and Anderson 1999) following the recommendation of Kierepka and Latch (2015) to combine statistical tests to assess IBD. Mantel tests were assessed using 5000 permutations in Arlequin (ver. 3.5.1.2, see http://cmpg.unibe.ch/software/arlequin35/, accessed 18 December 2013; Excoffier and Lischer 2010). Distance-based (db)RDA used the Fstat matrix against the sampling locations run through the R package vegan (ver. 2.4–2, J. Oksanen, F. G. Blanchet, R. Kindt, P. Legendre, P. R. Minchin, R. O’Hara, G. L. Simpson, P. Solymos, M. H. H. Stevens, and H. Wagner, see https://cran.r-project.org/web/packages/vegan/index.html, accessed 7 February 2018). IBD analyses were conducted on the complete dataset and on the sampling locations from Terrigal north (based on the K = 2 outcome from the structure and pairwise D_EST analyses).

A population genetic self-assignment test using a Rannala and Mountain Bayesian method with threshold 0.05 (Rannala and Mountain 1997) was conducted in GENECLASS2 (Piry et al. 2004) to determine the probability of correctly assigning an individual to a stock in the K = 2 stock model. The test was run both including and excluding the mixed Eden sampling location from the analysis.

The linkage disequilibrium method for estimating N_e (LDN_e; Waples and Do 2008) was applied using NeEstimator 2.01 (Do et al. 2014) to the samples from each collection area. For each set of samples, all loci and individuals were used and low frequency alleles were discarded if their observed frequency was below P_crit, the minimum allowed allele frequency. P_crit was chosen based on the sample size for the particular area according to standard methodology (Waples and Do 2010). Owing to the similar sample sizes from each area, P_crit = 0.02 for all the LDN_e estimates produced. For each estimate, a 95% confidence interval (CI) was also calculated according to the revised jack-knife method of Jones et al. (2016).

In all, 449 C. auratus were collected and genotyped as part of this study. Samples were obtained from nine collection areas over a 2278-km range from Rockhampton in the north to Tas. and Geelong in the south (Fig. 1). Sampled fish ranged in size from 150 to 810 mm fork length (FL), with an approximately equal sex ratio (sex recorded in ~50% of samples). Genotypes were obtained for all animals, across nine loci, with missing values per locus ranging from 0.2 to 9.4% (Table 1). Amplification difficulties due to poor-quality DNA samples were improved by re-extracting and eluting samples into a smaller volume (50 instead of 100 µL), then amplifying the difficult loci as singletons using undiluted DNA instead of multiplexed reactions. A relatedness screen of the samples identified a duplicate animal in the dataset. One of the duplicates was removed from subsequent analyses (taking the total number of C. auratus screened to 448; Table 1).

Table 1.  Summary statistics of nine microsatellite loci for Australian east coast Chrysophrys auratus
n scored, number of individuals genotyped; A, number of alleles; H_Obs, observed heterozygosity; H_Exp, expected heterozygosity; PIC, polymorphic information content

High allelic diversity was observed in the data, with 18 or more alleles found in two-thirds of the loci (allele number per locus ranged from 8 to 27; Table 1). Observed heterozygosity ranged between 0.342 and 0.879 and was lower than expected for all loci. With Bonferroni correction, only one locus by sampling location comparison was out of HWE (Sunshine Coast with Sal10; Table 2). The inbreeding coefficient (F_IS) value for this comparison was positive, indicating a heterozygote deficiency, most likely due to scoring errors. Although P-values were low for some of the other exact tests, there was no strong pattern of deviation from HWE by locus or by location; however, we note that deviations at one location (Sunshine Coast) and one locus (Pma68–23) may not be due to chance. Running the code (x <- replicate(1000000, sum(runif(9) < 0.05)) mean(x > = 3) mean(x > = 2)) in R determined the chance of obtaining three P-values <0.05 from a uniform distribution (0–1) if nine values were selected randomly a million times. The result was extremely unlikely due to chance, suggesting something biological is possibly happening at this site. To evaluate the possible effect of these deviations on population genetic structure, the structure analysis was repeated without Sunshine Coast and locus Pma68–23 (see Fig. S1). The results did not change the overall outcome, or interpretation, suggesting the genetic signal in the data was robust to the deviations. As a result, all loci and collection sites were retained for subsequent analyses.

Table 2.  Population × locus summary P-values for exact tests assessing Hardy–Weinberg equilibrium (HWE) with Bonferroni corrected α = 0.0056 for nine tests
The number of animals genotyped is given in parentheses (total n = 448). Significant P-values are in bold

Screening for LD using exact G-tests (log-likelihood ratio tests) to assess each pair of loci within each collection location identified three significant pairs (Sunshine Coast: Pma1 with CM003195; Eden: GA2A with Pma68–23; Eden: Pma22–9 with Pma68–23). The apparent genetic association of these pairs was not supported in other sampling locations and is likely an artefact of reduced sample numbers at these sites due to missing genotypes (Gordon et al. 1999; Akey et al. 2001). For this reason, at the conclusion of the exploratory analyses, all nine loci were determined to be suitable (largely meeting the assumptions of HWE and no LD) for downstream analyses. Based on the power analysis, the microsatellite data contain sufficient power to detect F_ST at or above 0.002 with 98% confidence (Fig. 2).

Fig. 2.  Power analysis results of simulations to assess the resolution of nine microsatellite loci to detect genetic structure in nine Chrysophrys auratus sampling locations given a range of divergence (fixation index, F_ST) levels. Probabilities reflect the average of 500 replicates for which the F_ST values were significantly different to zero using Fishers exact tests (POWSIM).

Analysing the data for genetic structure using Bayesian modelling (structure; Pritchard et al. 2000) identified a two-cluster model as the most likely fit for all models, with a genetic transition occurring around Eden on the NSW south coast (Fig. 3a). The optimal two-cluster model was determined by the highest average likelihood score (Fig. 3b) and the highest change in mean likelihood score (ΔK; Fig. 3c). The southern stock extends from Geelong, past Tas., then north to Lakes Entrance. The sampling locations flanking Eden (Terrigal, Eden, Lakes Entrance and Tas.) exhibit an intermediate level of mixing between the two stocks.

Fig. 3.  (a) Predicted genetic stock structure of Chrysophrys auratus based on population clustering of microsatellite data using a Bayesian model-based analysis. Vertical lines correspond to individual fish, which are coloured by the posterior probability proportions (log likelihood for each K, L(k)) of their genotype based on a K = 2 cluster model determined using (b) the best likelihood score at K = 2 and (c) the optimal K = 2 cluster model determined by the peak point on a plot of changes in mean likelihood scores (delta K) against cluster size (K). (a) Fish are plotted from north (left) to south (right) by sampling location. Rock, Rockhampton (Qld); Sun Cst, Sunshine Coast (Qld); Coffs H, Coffs Harbour (NSW); Wallis Lk, Wallis Lake (NSW); Lks Ent, Lakes Entrance (Vic.); Tas., Tasmania. The order of pooling for population pairwise comparisons in Table 3 are shown below the plot.

Using pairwise D_EST statistics to pool undifferentiated adjacent collection locations, the spatial boundaries of potential genetic stocks were further investigated (Table 3a). As a mixed transition zone between the northern and southern stocks, Eden was excluded from pooling. The first round of pooling grouped Rockhampton with Sunshine Coast, Coffs Harbour with Wallis Lake, and Lakes Entrance with Tas. (note, Lakes Entrance is a linear neighbour to both Tas. and Geelong). The second round of pooling combined Coffs Harbour + Wallis Lake with Terrigal in the north and Lakes Entrance + Tas. with Geelong in the south. The third round of pooling combined Rockhampton + Sunshine Coast with Coffs Harbour + Wallis Lake + Terrigal (diagram summarising pooling order below Fig. 3a Table 3a). The end result of three rounds of pooling was two significantly different stocks (Table 3b), the first extending Rockhampton to Terrigal and the second extending south from Lakes Entrance to Geelong, including Tas. with a significantly different D_EST = 0.048 (P = 0.001).

Table 3.  Pairwise Jost’s estimate of differentiation (D_EST; Jost 2008) below the diagonal and P-values above the diagonal, based on nine microsatellite loci from 448 individuals of Chrysophrys auratus among individual sampling locations (a) (significant Bonferroni-corrected P-values in bold, α = 0.0014) and pooled sampling locations (b) (significant Bonferroni-corrected P-values in bold, α = 0.017)
The final set of pooled locations was obtained after pooling strictly adjacent sampling locations that showed no significant pairwise D_EST until all pairwise D_EST were significantly different (as a transition zone, Eden was excluded from pooling; diagram at base of Fig. 3a summarises the order of pooling). Pooled sampling locations are as follows: Pool 1, Terrigal north, Rockhampton, Sunshine Coast, Coffs Harbour, Wallis Lake and Terrigal; Pool 2, Lakes Entrance south, Lakes Entrance, Tasmania and Geelong. Negative D_EST values were assigned zero. NSW, New South Wales; Qld, Queensland; Vic., Victoria

The optimum number of clusters (K with the lowest BIC score) determined using the DAPC find.clusters function was K = 1, although the BIC score for K = 2 was only a little higher (Fig. 4a). Results of the analyses run for a priori stock number K = 2, based on the outcomes from the structure analysis and population pairwise D_EST analysis, clearly showed a distinct northern and southern stock with a region of overlap at Eden (Fig. 4b).

Fig. 4.  Discriminant analysis of principal components (DAPC) using nine loci. (a) Plot of cluster size (K) v. Bayesian information criterion (BIC) to predict optimal K using scaling. (b) Percentage frequency plot of K = 2 stock a priori model separating north of Eden (northern cluster) from Eden south (southern cluster).

Genetic distance was found to be correlated with linear geographic distance using an IBD genetic model on the complete dataset (Fig. 5). The geographic distance between collection areas explained 34% of their genetic distance (R² = 0.3439). The IBD signal for the complete dataset was significant using both a Mantel test (5000 permutations; P = 0.0014) and the dbRDA analysis (P = 0.001). Focusing on the five sampling locations from Terrigal north, the northern stock in the K = 2 cluster model, IBD was no longer detected using either test (Mantel, P = 0.27; dbRDA, P = 0.217).

Fig. 5.  Correlation of genetic distance (calculated as D_EST/(1 – D_EST), where D_EST is Jost’s estimate of differentiation) and geographic coastal distance (isolation by distance (IBD) model) based on nine microsatellite loci and 448 Chrysophrys auratus collected from nine sampling locations along the east coast of Australia.

Loci Pma1, GA2A and CM000278 were excluded from the AMOVA because they were above the missing data threshold of 5%. Based on the structure and D_EST results, Eden was excluded from the analysis because of its mixed nature, and a K = 2 stock model was tested. Analysis of the remaining eight sampling locations and six loci using AMOVA found significant differentiation to support the two-stock model, although stock accounted for only 1.22% of the genetic variability (among group covariance F_CT = 0.01219, P < 0.0001). Most of the genetic variation in C. auratus (98.73%) was explained by within-sampling location differences.

Population genetic assignment tests correctly assigned 76.3% (including Eden) and 81.7% (excluding Eden) of samples to the correct stock in a K = 2 model. There was no obvious bias in the direction of incorrect assignments, with 18% of the northern stock incorrectly assigned to the south and 19% of animals from the southern stock incorrectly assigned to the north (K = 2 model excluding Eden).

Using the linkage disequilibrium method, estimates of LDN_e were calculated from six of the nine collection areas (Table 4). It was not possible to resolve LDN_e for Rockhampton, Wallis Lake or Geelong, possibly due to the fairly low power of the estimator and small sample sizes analysed. The upper 95% CI was infinite (i.e. could not be estimated) for any of the sites. Terrigal returned an LDN_e estimate considerably higher than any other site (7- to 40-fold larger). The wide range of sampling in Tas. may also be lowering the LDN_e estimate from this region.

Table 4.  Estimates of effective population size calculated using the linkage disequilibrium method (LDN_e), with upper and lower confidence intervals (CI) estimated using jack-knifing, for east coast Chrysophrys auratus sampling locations using a minimum allowed allele frequency (P_crit) of 0.02
Non-measurable values are indicated by a dash

Testing the null hypothesis of one continuous genetic stock of snapper along Australia’s eastern coast is important to ensure that the most appropriate spatial structure is applied to assessment and management. This study, using microsatellite markers, has confirmed a genetic subdivision of C. auratus populations along the east coast that was previously suggested based on a weak allozyme signal driven by a single locus (Sumpton et al. 2008). Temporal replication of genetic sampling spanning more than a decade has also offered a unique opportunity to investigate the geographic stability of this genetic disjunction.

The present study has revealed a consistent genetic signal, across all methods of statistical analysis (Bayesian modelling, DAPC and F_ST/D_EST statistics) for two distinct genetic stocks. The two genetic stocks transition around Eden in southern NSW. The northern stock spans most of NSW and Qld, and incorporates the region where most of the east coast catches are derived. The southern stock is most distinct in Geelong, but continues through Bass Strait, past Tas., extending up the east coast to Eden, with a small amount of encroachment north to Terrigal. The earlier study by Sumpton et al. (2008) also identified an east coast genetic disjunction but, in that study, it was identified over 400 km further north than in the present study, between Sydney and Forster on the central coast of NSW (Fig. 1). If the genetic disjunction is a consistent feature of this species on the east coast, then the genetic disjunction of Sumpton et al. (2008), based on samples collected in the mid-1990s, is likely the same genetic transition zone identified here, with the shift reflecting a southward movement of the ranges of the two genetic stocks.

The shifting genetic transition zone of the C. auratus stocks may be linked to ocean currents, temperature and possibly salinity. Long-term ocean temperature monitoring shows that the southward penetration of the EAC has increased over the past 60 years, resulting in a poleward advance of warmer and saltier water (Ridgway 2007). Water temperature has been shown to be linked to spawning periods and spawning success in C. auratus (Francis 1993; Lenanton et al. 2009). The southward-shifting EAC has been associated with long-term shifts in the abundance and distribution of other temperate fish species (Last et al. 2011).

The mechanism responsible for creating distinct northern and southern stocks may also be linked to water temperature. Fish living in warmer waters at the northern end of the range spawn during winter, whereas fish living in more temperate waters spawn ~3 months later in spring–summer (Ferrell and Sumpton 1997). As a result, northern fish have an extended growing season and they mature earlier than their southern counterparts (Ferrell and Sumpton 1997; Stewart et al. 2010). Biological parameters for the genetic stock occurring south of Eden are not well characterised; however, there are known spawning aggregations on inshore reefs off eastern Vic. in November, and recruitment of postlarvae and small juveniles occurs in eastern Victorian estuaries in late summer and autumn (Hamer and Jenkins 2004). The structure plot (Fig. 3a) shows that the two stocks overlap and transition around Eden. It is unclear whether the animals in this transition zone are interbreeding to form admixed animals (individual fish exhibit both northern and southern genotypes) or whether they co-occur in this region without interbreeding. There is a suggestion of admixture off Eden, as well as in neighbouring sampling locations, but this may be an artefact of the low number of microsatellite loci and level of missing data (this is the most likely explanation for the apparently admixed animals off the Sunshine Coast). Future studies should focus on screening additional genetic markers and sampling locations at a finer scale between Lakes Entrance and Terrigal to help to resolve this issue. Additional genetic markers may also reveal further genetic signals to differentiate Geelong from Lakes Entrance and Tas.

Standardised fixation measures may be less comparable between studies if mutation rates differ and heterozygosity is high (Leng and Zhang 2011). The total standardised fixation index D_EST (0.0232) estimate for the east coast samples was low but significant (P = 0.001). The estimate was almost an order of magnitude higher than the value recorded among Western Australian C. auratus populations from Shark Bay D_EST (0.002; Gardner et al. 2017), with the caveat being that only locus Pma1 overlapped between the two studies.

By pooling genetically undifferentiated adjacent collection areas, the spatial boundaries within the two east coast genetic stocks were further investigated using comparisons of pairwise D_EST. Following three rounds of pooling, two genetically distinct stocks were identified, the first extending from Rockhampton to Terrigal and the second spanning from Lakes Entrance to Geelong. Prior to pooling, Eden was the only sampling location with no significant D_EST values, a result consistent with it having a mixture of both the northern and southern stocks.

A priori expectations based on tagging data were that Geelong would separate from east coast C. auratus stocks (Sanders 1974; Hamer et al. 2011). The microsatellite data do not support this separation. This result is consistent with allozyme results that also failed to differentiate between Victorian sampling locations spanning from Portland in the west to Lake Tyers in the east (Meggs et al. 2003). Although not significant, population pairwise D_EST values within the southern stock were greater than those measured between sampling locations in the northern stock. Tasmanian samples were somewhat intermediate genetically between Geelong and Lakes Entrance. Low catch numbers unfortunately resulted in the Tasmanian samples being represented by fish collected from both north-west and eastern waters (see Fig. 1). The Tasmanian sample may have captured more than one genetic stock. With increased sampling, and additional genetic markers, it is possible that genetic structure in southern Australian C. auratus will become apparent.

High levels of connectivity were found among collection locations with a weak IBD signature detected using only the complete dataset. With increasing geographic distance, a linear increase in genetic differences was observed. The slope of the IBD correlation, 7 × 10^–6, falls within the interquartile range reported for fish stocks based on a meta-analysis of marine fish (Cooke et al. 2016). The slightly positive slope value indicates that fish are probably not actively recruiting back to natal sites (Cooke et al. 2016), but neither are they moving enormous distances; they are sharing their DNA with neighbouring locations, either by mixed spawning aggregations or larval dispersal. The IBD signature could not be detected when the microsatellite dataset was reduced to the northern stock alone (five sites from Terrigal north) although a weak signal was found by Sumpton et al. (2008) using allozyme analysis. Biological knowledge of the fish would suggest that an IBD pattern likely does exist in the northern stock but has not been detected with the current genetic markers. A similar result was obtained using a microsatellite analyses of New Zealand C. auratus over a 900-km range (Ashton 2013), with stocks in that study found to be largely panmictic with no IBD signal and with low-level genetic differentiation between sites.

This study has demonstrated that suitable genetic markers exist to characterise the population genetic structure of C. auratus stocks in Australian waters. Further research is needed to characterise stocks in southern and western Australian waters. A broad-scale nationwide study is recommended, with a consistent set of variable genetic markers, before fine-scale spatial variation is assessed to assist in the management of local stocks. It would also be valuable to compare Australian C. auratus stocks to populations from New Zealand, which were found to be genetically distinct using allozymes (Meggs et al. 2003). Although the microsatellite loci screened in this study were the same as those used by Le Port et al. (2017), unfortunately different size-scaling ladders were used by the two laboratories for scoring alleles, and consequently the results require cross-validation before they can be combined.

It would be premature to infer population numbers (N) from the LDN_e estimates calculated without a better understanding of how the two numbers correlate in C. auratus, but some interesting observations can be made from the results. It is promising that LDN_e values were estimated using a dataset limited by low sample numbers, few loci and with a high occurrence of rare alleles (which are currently excluded by the P_crit). Having infinite upper bounds, the estimates are currently fairly meaningless; however, with more intensive sampling and the addition of more loci, the LDN_e estimates could offer a comparative means to assess trends in the abundance of local stocks.

Microsatellite analysis supports a two-stock genetic model for C. auratus along Australia’s east coast. The northern stock extends from Rockhampton to Eden; the northern and southern stocks overlap around Eden and the southern stock then extends south from Eden, past Lakes Entrance and Tas., extending westward to Geelong. The southern stock will likely be found to extend further west. Other studies have suggested that mixing of snapper between stocks to the east and west of Wilsons Promontory in Vic. is limited (Coutin et al. 2003; Hamer et al. 2011). The present study has not detected a significant barrier to gene flow to separate Geelong (Port Phillip Bay) samples from eastern Victorian or Tasmanian samples. The consequence of this finding, in terms of stock assessment and fisheries management along the east coast, is that Qld and NSW fisheries should be assessed together. Preliminary genetic analyses suggest that, with increased sampling, the size of local C. auratus stocks may be able to be independently estimated using DNA tissue samples. This finding could prove valuable for the local management of populations.

Jennifer Ovenden is an Associate Editor for Marine and Freshwater Research. Despite this relationship, she did not at any stage have Associate Editor-level access to this manuscript while in peer review, as is the standard practice when handling manuscripts submitted by an editor to this Journal. Marine and Freshwater Research encourages its editors to publish in the journal and they are kept totally separate from the decision-making processes for their manuscripts. The remaining authors have no conflicts of interest to declare.

This project was supported by the Australian Fisheries Research and Development Corporation (Project number 2015/216).