Genetic structure and recent population expansion in the commercially harvested deep-sea decapod, Metanephrops challengeri (Crustacea: Decapoda)

ABSTRACT

Determining the genetic structure and population dynamics of harvested species are key aspects of effective fisheries management. New Zealand scampi (Metanephrops challengeri) is a species of nethropid lobster that supports a commercial fishery on the continental shelf and slope around New Zealand. To determine the genetic structure and historical population demography of New Zealand scampi, we sequenced a 623 base-pair fragment of the mitochondrial gene cytochrome c oxidase subunit one (CO1) from 321 individuals across six widely separated sites off the coast of mainland New Zealand and the subantarctic Auckland Islands, representing much of the species’ distribution. The Auckland Islands population was found to be genetically distinct from all other sampled populations, potentially resulting from an absence of continuous suitable habitat between the central-eastern South Island and the Auckland Islands. An isolation-by-distance pattern of genetic structure was also detected. With the exception of the Mercury Islands, all sampled populations showed signs of recent population expansion, potentially linked to changes in habitat availability during the glacial-interglacial cycles of the Pleistocene. We recommend that the current NZ scampi management areas and harvesting rates for each remain unchanged.

KEYWORDS:

• Scampi
• Nephropidae
• fisheries
• phylogeography
• stock
• New Zealand

Introduction

Identification of fisheries stock units is an essential aspect of effective fisheries management. The overall goal of fisheries management is to balance harvest-induced mortality with natural recruitment and immigration to maintain a self-sustaining population (Reiss et al. 2009). Small populations, or those that have undergone a bottleneck, face increased risk of extirpation and can hold a lower level of genetic diversity, reducing the adaptive potential of the population. The effect of genetic drift increases with decreasing population size, which reduces genetic diversity, and, in some cases, causes deleterious alleles to increase in frequency (Laikre et al. 2005). Dispersal and interbreeding between populations can influence genetic connectivity and demographic stability (Lowe and Allendorf 2010; Faurby and Barber 2012). The delineation of stock boundaries within the marine environment can be difficult. Patterns of population connectivity in deep-sea species remain poorly understood (Hilário et al. 2015; Baco et al. 2016), and the majority of marine phylogeographic studies to date primarily consider coastal species due to the difficulty and cost of effectively sampling the deep-sea (Bors et al. 2012). Deep-sea species generally have higher levels of population connectivity compared to shallow-water species (Baco et al. 2016), due, in part, to the relatively few barriers to dispersal in the deep-sea environment (Baco et al. 2016).

For many marine invertebrates, dispersal primarily occurs during their pelagic larval phase (Cowen and Sponaugle 2009), the duration of which (pelagic larval duration; PLD) is often used as a proxy for dispersal potential (Marko and Hart 2017). Direct study of larval dispersal is challenging because of the vast areas and large numbers of minute larvae involved (Levin 2006; Cowen and Sponaugle 2009). Moreover, larvae may not persist for long periods in the pelagic environment. This has meant that indirect methods (e.g. analyses of genetic structure) are frequently used to estimate the extent of larval dispersal. Complex oceanographic conditions and variation in larval behaviour can also influence the distance and success rate of dispersal. Ocean currents can enable dispersal far away from natal areas (Cowen et al. 2000), or, if larvae become trapped within local eddies, promote self-recruitment (Lobel and Robinson 1986). Halo- and thermo-clines can also restrict movement of marine organisms or form semi-permeable barriers (e.g. Jørgensen et al. 2005; Thornhill et al. 2008). Variation in the distributions and connectivity patterns among and within species is due to differing habitat requirements and availability, physiological tolerances to a range of abiotic factors, life-history strategies, and chance events (Ross et al. 2009). With so many biotic and abiotic factors involved, it is often misleading to apply generalisations about population connectivity from one species to another (e.g. Goldstien et al. 2006; Fleming et al. 2018). Therefore, patterns of population connectivity need to be assessed separately for each species (Wennerström et al. 2013).

New Zealand scampi (Metanephrops challengeri; NZ scampi hereafter) are targeted by a commercial trawl-based fishery, with the primary fishery areas comprising the region east of Great Barrier Island to the western Bay of Plenty in the northern North Island, east of the Wairarapa in the southern North Island, the Chatham Rise east of mainland New Zealand, and off the Auckland Islands, 450 km south of the mainland (Ministry for Primary Industries 2016). The fishery is currently managed as eleven separate stocks within the New Zealand Exclusive Economic Zone (EEZ) under New Zealand's quota management system. NZ scampi are endemic to New Zealand and inhabit soft benthic substrates on the continental shelf and slope at depths of 140–670 m (Yaldwyn and Webber 2011), in highest abundance between 200 and 500 m (Smith 1999). NZ scampi have a bipartite life cycle comprising a short pelagic larval stage of 3–5 days (Wear 1976) before benthic settlement. A mark-recapture study of adult individuals supported the current stock structure model by suggesting that adults do not traverse large distances, possibly moving as little as 500 m post settlement (Cryer and Stotter 1999). The short pelagic phase and limited movement of adults suggest that NZ scampi have low dispersal ability and predicted low levels of gene flow amongst geographically separated populations. However, there has been limited research into the genetic structure of NZ scampi populations, which would enable the efficacy of the current proposed stock structure to be tested.

Low-resolution allozyme data reported that the Auckland Islands population was genetically distinct from all other regions on the eastern mainland, with significant differences in allele frequencies at three of five polymorphic loci, and that there was significant genetic differentiation between the Bay of Plenty and the Chatham Rise populations (Smith 1999). However, this genetic differentiation was weaker, with only one polymorphic loci exhibiting significantly different allele frequencies between the latter populations (Smith 1999). The aim of the present study was to use DNA sequences from the faster evolving mitochondrial gene cytochrome c oxidase subunit 1 (CO1) to test for phylogeographic structure, and describe genetic diversity and past demography among and within NZ scampi populations. These data could then inform management of the NZ scampi fishery.

Methods

Sample collection

Samples (Table 1) were collected via trawling in Bay of Plenty, south of Cape Kidnappers and off Wairarapa during a National Institute of Water and Atmospheric Research (NIWA) NZ scampi survey (KAH1501). Samples collected from off the Auckland Islands were taken during a separate survey (KAH1603), and Chatham Rise samples were sourced from commercial fishing vessels (see Figure 1 for a map of sampling localities). Trawls were conducted between 300 and 500 m water depth. A section of muscle tissue was excised from the tail and stored at 4°C in 80% ethanol.

Figure 1. Map of New Zealand indicating sampling locations and areas of interest. AI: Auckland Islands; CR: Chatham Rise; CP: Castlepoint; CK: Cape Kidnappers; WI: White Island; MI: Mercury Islands. Latitude, longitude, 200 m isobath, and 1000 m isobath shown.

Table 1. Sampling information, summary statistics, estimates of genetic diversity, and neutrality test statistics for 321 partial NZ scampi CO1 sequences. For neutrality test values (Tajima's D and Fu's F-statistic), bold values indicate statistical significance. Key: n: Sample Size, S: Segregating Sites, T_S: Transitions, T_V: Transversions, H_N: Number of Haplotypes, H_P: Number of Private Haplotypes, H_D: Haplotype Diversity, π: Nucleotide Diversity (standard deviation), k: Average Number of Nucleotide Substitutions.

Sample site	Lat, long	n	S	Ts	Tv	HN	HP	HD	π	k	Tajima's D (p-value)	Fu's F-statistic (p-value)
Mercury Islands	36°20′S, 176°12′E	63	41	29	12	25	16	0.882	0.00724 (0.0008)	4.479	−1.604 ( 0.027)	−9.819 ( 0.004)
White Island	37°36′S, 177°10′E	48	30	24	7	24	15	0.922	0.00629 (0.0005)	3.889	−1.421 (0.058)	−13.123 ( <0.000)
Cape Kidnappers	39°52′S, 177°20′E	63	23	17	7	19	9	0.839	0.00541 (0.0006)	3.348	−0.984 (0.168)	−6.145 ( 0.017)
Castlepoint	40°53′S, 176°39′E	55	30	19	11	32	21	0.919	0.00680 (0.0006)	4.215	−1.212 (0.102)	−19.343 ( <0.000)
Chatham Rise	43°24′S, 176°03′E	46	24	18	6	21	10	0.850	0.00581 (0.0005)	3.603	−1.05 (0.148)	−8.775 ( 0.001)
Auckland Islands	51° 12′S, 166°41′E	46	20	14	7	22	21	0.885	0.00345 (0.0004)	2.136	−1.715 ( 0.022)	−18.649 ( <0.000)
All Sites		321	79	55	27	107	92	0.914	0.00810 (0.0003)	5.004	−1.766 ( 0.007)	−24.952 ( <0.000)

DNA extraction, polymerase chain reaction, and DNA sequencing

DNA was extracted using a standard phenol–chloroform method (Sambrook et al. 1989), quantified using a Nanodrop ND1000 spectrophotometer (Thermo Fisher Scientific), and diluted to 50 ng/µL in Tris-EDTA buffer (10 mM Tris-HCl pH 8.0, 1 mM EDTA). A 710 base-pair (bp) fragment of CO1 was amplified using the polymerase chain reaction (PCR) with the universal invertebrate primer pair LCO1490 (5′-GGTCAACAAATCATAAAGATATTGG-3′)/HCO2198 (5′-TAAACTTCAGGGTGACCAAAAAATCA-3′) (Folmer et al. 1994). Each PCR (25 µL) contained 1 µl of DNA, 1x Bioline PCR Buffer (67 mM Tris-HCl, 16 mM (NH₄)₂SO₄), 3 mM MgCl₂, 0.06 mM Bovine Serum Albumin, 0.1 µM each of forward and reverse primer, 0.4 mM dNTPs, and 1 unit of Biotaq DNA Polymerase (Bioline Reagents Ltd). Thermocycling was carried out on a Biometra TProfessional Thermocycler (Göttingen, Germany) using the following parameters: 95°C 3 min, 40 cycles of 95°C 35 s, 50°C 35 s, and 72°C 45 s, with a final extension of 72°C 10 min. PCRs were visualised on a 2% agarose gel stained with ethidium bromide and exposed to UV light. Successful reactions were then purified using ExoSAP-IT (Amersham Parmacia Biotech) following the manufacturer's instructions. DNA sequencing was conducted by the Macrogen Sequencing Service (Macrogen Inc. Seoul, South Korea) using a 3730xl Genetic Analyzer (Applied Biosystems). DNA sequences were checked by eye, edited where appropriate, trimmed to 623 bp and aligned with Geneious (Version 8.0.5) (Kearse et al. 2012) using default parameters.

Genetic diversity and population structure

Summary statistics and estimates of genetic diversity were calculated using DnaSP (Version 5.10) (Librado and Rozas 2009) and Arlequin (Version 3.5.2.2) (Excoffier and Lischer 2010). Summary statistics included the number of segregating sites (S), frequency of transitions (T_s) and transversions (T_v), nucleotide diversity (π), haplotype diversity (H_D), and the average number of nucleotide differences between sequences (k). Population structure was assessed using Bayesian Analysis of Population Structure (BAPS) (Version 6.0), with the implementation of the linked loci clustering model (Corander and Tang 2007) and multiple upper bounds to the number of clusters (k). K values were set from 2 to 10 with each tested 5 times to determine the most likely number of clusters. PopART (Leigh and Bryant 2015) was used to generate a median joining haplotype network (Bandelt et al. 1999), to enable haplotype frequencies and relationships to be visualised.

Pairwise ɸ_ST values and AMOVAs (Analysis of Molecular Variance) were calculated within Arlequin using haplotype frequencies. Each analysis was bootstrapped with 10,000 permutations and a Bonferroni correction applied to pairwise ɸ_ST estimates to reduce the chance of type one errors (false positives). AMOVAs were used to test support for different scenarios of population genetic structure through the a priori assignment of DNA sequences from different populations into groups. Three scenarios of population structure were constructed and tested using an AMOVA: (1) DNA sequences from different populations were placed into two groups, one consisting of the Auckland Islands, the other comprising all other populations. This grouping was chosen based upon the results of the BAPS analysis and haplotype network; (2) Populations were placed into three groups, one comprising the four North Island sites, the second group comprising only the Chatham Rise and the third, the Auckland Islands; (3) Populations were placed into two groups, one comprising all sites north of 42°S, the other comprising the Chatham Rise and Auckland Islands. This is based on previously published reports of genetic disjunctions in marine invertebrates around this latitude (e.g. Waters and Roy 2004; Ayers and Waters 2005; Goldstien et al. 2006; Veale and Lavery 2011).

The Isolation-By-Distance Web Service (http://ibdws.sdsu.edu/~ibdws/) (Jensen et al. 2005) was used to test for a correlation between genetic differentiation and geographic distance (Isolation-by-Distance). A Mantel's test with 1,000 randomisations was used on the regression analysis to assess significance. Minimum over-water distances (km) between sampling locations were estimated using Google Earth (version 7.1). Genetic distances were transformed using Slatkin's linearised genetic differentiation formula [ɸ_ST / (1−ɸ_ST)] (Slatkin 1995), following Rousset (1997), who suggested using Slatkin's linearised distance and geographic distance for one-dimensional habitats. A one-dimensional habitat was defined as habitat where the distance between sampled populations is greater than the length of the habitat (Rousset 1997). Habitat within this study was defined as one-dimensional due to the large latitudinal range sampled. A final regression analysis was conducted that excluded the Auckland Islands to check if the relatively high geographic- and genetic distances between the Auckland Islands and other locations strongly influenced the previous analyses.

Demographic history

Multiple approaches were used to examine the demographic history of NZ scampi. These include calculation of Tajima's D (Tajima 1989), and Fu's F statistics (Fu 1997), mismatch distributions of pairwise DNA sequence divergence, and Bayesian skyline plots (Drummond et al. 2005). Tajima's D and Fu's F-statistics were calculated in Arlequin and bootstrapped with 10,000 permutations. Tajima's D was deemed to be statistically significant at p-values of <0.05, and Fu's F-statistic was deemed statistically significant at p-values of <0.02 (Fu 1997). Mismatch distributions (Rogers and Harpending 1992) were calculated for each separate population using DnaSP. Bayesian skyline plots were constructed in BEAST v1.8.2 (Drummond et al. 2012). Initially, two separate analyses were run; one containing the Auckland Islands sequence data, and another containing all other DNA sequences, to meet the Bayesian skyline plot assumption of panmixia. However, preliminary results suggested that the DNA sequence data of the Auckland Islands population alone did not contain enough power to reject the assumption of constant population size through time (data not shown), therefore, further Bayesian skyline plot analysis of the Auckland Islands population was not attempted. For the remaining DNA sequences (those obtained from populations along the eastern coast of the North Island and the Chatham Rise) multiple analyses were run using a strict clock model with the following mutation rates: 1 × 10⁻⁸, 7 × 10⁻⁹, and 1.3 × 10⁻⁸ substitutions/site/year. These mutation rates correspond to an intermediate estimate of decapod CO1 mutation rate (Zacarias 2013), alongside a rates reported for snapping shrimp within the genus Alpheus (Knowlton and Weigt 1998), and Jamaican land crabs (Schubart et al. 1998), respectively. The TN93+I+G model of nucleotide substitution was used as it was identified as the most appropriate model using the AIC criterion in jModelTest v2.1.1 (Darriba et al. 2012). A total of 50⁶ iterations were analysed and the first 10% discarded as burn-in. The resulting trace file was viewed in Tracer v1.7 (Rambaut et al. 2018), and the effective sample sizes of priors and posterior probabilities were confirmed to be greater than 200.

Results

Genetic diversity and population structure

The final sequence alignment comprised 321 sequences with a total of 107 haplotypes. There were 79 variable sites, of which, 39 were parsimony informative and 40 were singletons. All populations contained private haplotypes, with the highest number of private haplotypes being detected from the Auckland Islands and Castlepoint sampling sites (H_P = 21 for both sites). All populations had low levels of nucleotide diversity (π = 0.00345–0.00724) (Table 1), and high levels of haplotype diversity (H_D = 0.839–0.922). The Mercury Islands site had the highest level of nucleotide diversity and highest average number of nucleotide substitutions (π = 0.00724, k = 4.479). The Auckland Islands samples featured almost exclusively private haplotypes (H_N = 22, H_P = 21), and contained lower levels of nucleotide diversity (π = 0.00345) and a lower number of nucleotide substitutions (k = 2.136) than the other sample sites (Table 1).

BAPS analysis indicated that there were five genetic clusters among the six sampled populations. Samples from the Auckland Islands formed one group, with a few samples from Cape Kidnappers (eastern North Island) and White Island (northern North Island) clustering with them (Figure 2). The remaining samples from the populations other than the Auckland Islands were intermixed among the other four clusters. The median joining haplotype network also recovered the same population structure as BAPS from the dataset. The haplotype network (Figure 3) exhibited four expanded starburst-like haplogroups, indicative of recent demographic expansion. All samples from the Auckland Islands belonged to a single haplogroup, along with some samples from Cape Kidnappers and Mercury Islands (Figure 3), separated by at least five substitutions from the other major haplogroups.

Figure 2. Bayesian analysis of population structure with k = 5. A total of 321 samples are represented. Colours represent different genetic clusters defined by BAPS. MI = Mercury Islands; WI: White Island; CK: Cape Kidnappers; CP: Castlepoint; CR: Chatham Rise; AI: Auckland Islands.

Figure 3. Median-joining haplotype network of 321 partial NZ scampi CO1 sequences. Each hatch mark represents one mutational step, while small black circles represent hypothetical unsampled haplotypes.

AMOVA analyses failed to detect a statistically significant difference among groups (ɸ_CT) in all tested scenarios (p = 0.13129–0.32337) (Table 2). However, statistically significant differences were detected among all populations (ɸ_ST) (p < 0.001) and within predefined groups (ɸ_SC) (p < 0.001) for all sample groupings. Most genetic variation within the dataset could be explained by the variance among populations (88.18%–94.67%), with small amounts explained by the differences within and among groups. Pairwise ɸ_ST values show statistically significant genetic differentiation between samples from the Auckland Islands and all other populations (Table 3). Small, yet statistically significant pairwise ɸ_ST values were also obtained when comparing sampling localities within Bay of Plenty (Mercury Islands and White Island) to those along the eastern coast of the North Island, and the Chatham Rise (Table 3). Of these, only the comparisons between the Mercury Islands and Cape Kidnappers, and the Mercury Islands and Chatham Rise remained statistically significant after Bonferroni correction.

Table 2. Analysis of Molecular Variance (AMOVA) analyses of M. challengeri using 321 CO1 sequences with 3 a priori geographic groupings. (1) AI, MI-CR: Group 1: Auckland Islands, Group 2: Chatham Rise, Castlepoint, Cape Kidnappers, White Island, Mercury Islands. (2) AI, CR, MI-CP: Group 1: Auckland Islands, Group 2: Chatham Rise, Group 3: Castlepoint, Cape Kidnappers, White Island, Mercury Islands. (3) AI+CR, MI-CP: Group 1: Auckland Islands, Chatham Rise, Group 2: Castlepoint, Cape Kidnappers, White Island, Mercury Islands. d.f. = degrees of freedom, Var = variance component, %Var = variance component in percentage of total variation.

Groups	Among groups					Among populations within groups					Within populations
	d.f	Var	%Var	ɸCT	P-value	d.f	Var	%Var	ɸSC	P-value	d.f	Var	%Var	ɸST	P-value
AI, MI-CR	1	0.0502	10.12	0.1012	0.16901	4	0.00847	1.71	0.01899	0.00069	315	0.43751	88.18	0.11824	<0.00000
AI, CR, MI-CP	2	0.02622	5.55	0.05545	0.13129	3	0.00910	1.93	0.02038	0.00069	315	0.43751	92.53	0.07470	<0.00000
AI+CR, MI-CP	1	0.00260	0.56	0.00563	0.32337	4	0.02203	4.77	0.04795	<0.00000	315	0.43751	94.67	0.05331	<0.00000

Table 3. Pairwise ɸ_ST values determined by CO1 haplotype frequencies. Pairwise ɸ_ST values were determined using 321 CO1 sequences and the haplotype frequencies between sampling localities. These are situated below the diagonal while the p-value associated with each comparison is located above the diagonal. Bold values indicate statistical significance after Bonferroni correction.

	Auckland Islands	Chatham Rise	Castlepoint	Cape Kidnappers	White Island	Mercury Islands
Auckland Islands	–	<0 . 0000	<0 . 0000	<0 . 0000	<0 . 0000	<0 . 0000
Chatham Rise	0 . 14058	–	0.70003	0.55994	0.01940	0 . 00158
Castlepoint	0 . 11063	0.0000	–	0.51866	0.02861	0.00554
Cape Kidnappers	0 . 13850	0.0000	0.00208	–	0.06009	0 . 00149
White Island	0 . 10273	0.04614	0.01786	0.01514	–	0.21127
Mercury Islands	0 . 09644	0 . 03791	0.03946	0 . 03878	0.00480	–

Mantel tests using Slatkin's linearised distance and geographic distance show a strong positive correlation between genetic and geographic distance (Figure S1) (R = 0.81–0.87). Interestingly, this statistically significant correlation increased slightly upon the removal of the Auckland Islands samples from the dataset (Figure S1B) (R = 0.8724; p = 0.01). In contrast, when the entire dataset was considered, a slightly weaker, non-statistically significant correlation was detected (R = 0.81, p = 0.052, Figure S1A). Large correlation coefficients indicate that geographic distance may explain approximately 80% of the variation in genetic distance between populations.

Demographic history

All neutrality tests returned negative values, though not all of these were statistically significant (Table 1). Fu's F-statistic was statistically significant for each population when analysed separately, and when all populations were combined. Tajima's D was statistically significant when considering all populations and when the Mercury Islands, White Island, and Auckland Islands were considered separately. A negative value from such neutrality tests can be due to an increase in population size or a bout of purifying selection (Tajima 1989). Mismatch distributions were consistent with the results returned by the neutrality tests and suggest that NZ scampi may have recently undergone a population expansion as some populations displayed genetic signatures of this event (Figure S2). Histograms for the Wairarapa, Cape Kidnappers, White Island and Mercury Islands sites, of the number of pairwise differences between sites, had a multimodal distribution. Conversely, histograms for the Auckland Islands and Chatham Rise sites had a unimodal distribution, consistent with recent population expansion. The Bayesian skyline plots estimated a stable effective population size over the past 500,000 years for all sampled populations from mainland New Zealand and the Chatham Rise (Figure 4).

Figure 4. Bayesian skyline plots with differing mutation rates reconstructing changes in NZ scampi effective population size through time. A: 7 × 10⁻⁹ substitutions/site/year. B: 1 × 10⁻⁸ substitutions/site/year. C: 1.3 × 10⁻⁸ substitutions/site/year. The blue-shaded areas represent 95% confidence intervals.

Discussion

Statistically significant levels of genetic differentiation were found between NZ scampi from the Auckland Islands compared to those from the Chatham Rise and eastern North Island combined. NZ scampi from these latter regions also exhibited small but significant levels of genetic differentiation. Similar results were reported by Smith (1999) based on allozyme data. The level of genetic structure between the Auckland Islands population and populations associated with mainland New Zealand is greater than reported levels of mitochondrial genetic structure between populations of European scampi (Nephrops norvegicus; Stamatis et al. 2004), similar to the genetic structure between populations of red spiny lobster (Panulirus interruptus; Perez-Enriquez et al. 2001), and significantly less than the genetic structure reported between populations of European spiny lobster (Palinurus elephas; Palero et al. 2008). The precise location of the genetic break between the Auckland Islands and eastern New Zealand reported here remains to be established due to the large distances between sample sites. In the vicinity of 45°S on the east coast of the South Island (along the Otago Shelf, Figure 1), there appears to be little co-occurrence of suitable habitat and bathymetry; potentially representing an impediment to gene flow. Our analyses suggest that NZ scampi populations have undergone a population expansion, although we cannot discern the specific timing of this event. This expansion may have been caused by the glacial-interglacial cycles of the Pleistocene, which similarly impacted lobster populations worldwide (see below).

Population structure

NZ scampi have a wide yet patchy distribution encompassing much of the continental shelf and slope off New Zealand (Smith 1999). Very few voucher specimens are held at the Museum of New Zealand Te Papa Tongarewa (R. Webber, pers. comm. 2017) or National Institute of Water Atmospheric Research (S. Mills, pers. comm. 2017). No voucher specimens are known from the east coast of the South Island between Timaru and Nugget Point, or east of Stewart Island, and neither are there records of commercial trawls targeting NZ scampi within these regions (Tuck 2014: Figure 1). There are no known records of NZ scampi off Otago Peninsula despite extensive trawling there (K. Probert, pers. comm. 2017). Off the eastern South Island, the continental shelf is relatively shallow, with only a few areas close to the upper bathymetric limit for NZ scampi (Mackay et al. 2015), and the continental slope drops off rapidly into unfavourable depths (Carter and Carter 1988). Land-derived gravels and sands make up much of the sediments within this area (Carter 1985), which may be unsuitable for burrow formation by NZ scampi. Suitable habitat for Metanephrops species is often correlated with sediment type and grain size (Wallner and Phillips 1990). The Otago region would have changed significantly during the Last Glacial Maximum, when eustatic sea levels were 120–130 m lower (Fleming et al. 1998) than at present. However, there would still have been little soft-sediment habitat within the bathymetric range of NZ scampi in this region during glacial maxima. Furthermore, the Clutha River on the east coast of the South Island has a strong influence on local ocean conditions such as temperature, salinity, and sediment load (Carter and Carter 1986; Hawke 1989), and the continental shelf and slope off Otago Peninsula is subject to density flows (Hawke 1989); either of which might present an impassable barrier to dispersing NZ scampi larvae or result in unfavourable habitat for adults.

Additionally, NZ scampi are not known to inhabit the waters off- or within- Fiordland (Brewin et al. 2008), possibly reflecting the unique hydrology of the region (Brewin et al. 2011). It seems probable that NZ scampi are either absent or in low abundance between 44°S and 46°20′S, and 45°30′S and 47°30′S, on the west and east coasts respectively of the South Island. Distribution disjunction could explain the marked genetic structure between the Auckland Islands and the other sampled populations. Additional sampling (e.g. the Challenger Plateau, west of Stewart Island, and the Antipodes Islands) could improve understanding of how and why NZ scampi populations are genetically structured.

Mantel tests showed a strong signal for genetic isolation by distance, although this relationship was not always statistically significant (Figure S1). The Auckland Islands are separated from the Chatham Rise by c. 1100 km, and there is significant genetic differentiation between NZ scampi populations in these areas. However, over a similar distance between the Chatham Rise and Mercury Islands, there was significantly less genetic structuring. While hypothesised phylogeographic barriers effectively isolate the Auckland Islands population, the near continuous population along the eastern coast of the North Island and the Chatham Rise may promote gene flow between neighbouring areas. Simulated larval dispersal models examining the subantarctic islands of New Zealand (including the Auckland Islands) demonstrate that larvae with a short PLD can only travel a short distance from their spawning site (Chiswell 2009), potentially limiting dispersal and gene flow between geographically separated areas.

Ocean currents are likely to impact the genetic and demographic connectivity of NZ scampi populations because their pelagic larvae are passive dispersers. However, the behaviour and positioning of larvae in the water-column are as yet poorly known, as is the effect and directionality of deeper currents in the bathymetric range of settled individuals. The main current systems that could influence the dispersal of individuals between NZ scampi populations examined in this study are the Tasman Front, which influences western New Zealand, especially in northern regions, and the Subtropical Front (STF), which extends to 500 m below sea level along the Otago shelf (Shaw et al. 1999), and passes to the south of the South Island before running eastwards along the Chatham Rise (see Chiswell et al. 2015). The STF has been proposed as a potential dispersal barrier to NZ scampi (Smith 1999), as populations on either side of the STF are genetically distinct (Smith 1999; this study). However, as above, the locality of the genetic disjunction cannot be precisely discerned due to the large geographic distance separating genetically distinct sampling sites. We observe no evidence for a north–south phylogeographic break around Cook Strait (Tables 2 & 3), as noted in studies of other marine invertebrates (e.g. Waters and Roy 2004; Ayers and Waters 2005; Veale and Lavery 2011). Many of the species exhibiting phylogeographic breaks in this area are coastal taxa, while the location of these breaks are often poorly defined and vary from species to species (Ross et al. 2009; Walton et al. 2019). Deep-sea species (i.e. NZ scampi) may be unconstrained by the biotic and abiotic factors that impact the connectivity of coastal species within this region.

Statistically significant genetic structure was detected amongst NZ scampi populations from the Bay of Plenty (Mercury Islands), the eastern coast of the North Island (Cape Kidnappers), and the Chatham Rise (Table 3). The East Auckland Current flows south along the north eastern coast of New Zealand, reaching the East Cape before much of this water is transported offshore. Offshore transport of water may influence the connectivity between populations of pelagic dispersers north and south of the East Cape by transporting pelagic larvae offshore and into waters that are unsuitable for larval settlement due to their depth, with populations of coastal marine invertebrates either side of the East Cape recorded as genetically distinct (Stevens and Hogg 2004; Veale and Lavery 2012). The current systems within these areas are further complicated by the presence of the East Cape and Wairarapa eddies, which have the potential to retain small drifting objects for up to 3 years (Chiswell and Roemmich 1998). These systems have been cited as a driver of self-recruitment in populations of Crustacea with an extremely long PLD (Chiswell and Booth 1999), but for species with a short PLD, such as NZ scampi, retention within these eddies will be more significant as mature larvae will be unable to reach habitat suitable for settlement (Ross et al. 2009). Populations of NZ scampi north and south of the East Cape are genetically connected, but genetic connectivity between populations does not imply that they are also demographically connected (Ovenden 2013). The East Cape could constitute a semi-permeable barrier for NZ scampi with successful dispersal occurring around it at low frequencies, leading to low levels of structure between populations either side of the East Cape.

Evolutionary insights determined via the sole analysis of mitochondrial DNA are limited by its maternal inheritance and lack of recombination. While our mitochondrial DNA results are supported by nuclear allozyme data (see Smith 1999), genomic analyses such as genotyping-by-sequencing (Elshire et al. 2011) may be able to detect additional genetic structure amongst NZ scampi populations, and further elucidate the recent evolutionary history of NZ scampi.

Demographic history

NZ scampi exhibit high levels of mitochondrial haplotype diversity. A total of 107 different CO1 haplotypes were identified from just 321 samples, which contained low levels of nucleotide diversity overall. These values are comparable to those exhibited by the CO1 region of M. mozambicus (Zacarias 2013), and slightly lower than observed levels in the deep-sea squat lobster Munida gracilis (Bors et al. 2012); the latter are syntopic with NZ scampi in parts of their respective distributions. High levels of haplotype diversity paired with low levels of nucleotide diversity suggests a demographic history of population size contraction before a recent expansion (Grant and Bowen 1998). Populations of triplefin fishes from southern New Zealand also showed high levels of haplotype diversity paired with low nucleotide diversity (Hickey et al. 2009). Additionally, a decrease in genetic diversity with increasing latitude has also been demonstrated in other species of crustacea (McMillen-Jackson and Bert 2003). The high effective population size of NZ scampi may reflect the relatively large area of available habitat, as greater habitat area is available to deep-water species than species constrained to coastal habitat. Mismatch distributions and the Bayesian skyline plot also reported patterns of population expansion.

The Pleistocene was characterised by a series of glacial-interglacial cycles that have left their mark upon the genetic population structure and demographic history of many taxa (see Hewitt 2000, 2004; Petit et al. 2003; Hickey et al. 2009). Pleistocene glacial cycles led to repeated extirpation-recolonisation events within the Southern Ocean fauna (Fraser et al. 2009; Nikula et al. 2010), and the onset of the present interglacial has allowed many Southern Hemisphere species to expand their distributions southward (Fraser et al. 2012). The northern and southern distribution limits of NZ scampi are bound by unfavourable bathymetry such that it is difficult to determine its thermal tolerance range. Southward isotherm tracking following the LGM, as recorded for a variety of marine invertebrates from shelf and slope depths (e.g. Maurea blacki, Marshall 1995; Zygochlamys delicatula, Beu 1999), could explain the genetic signal of population expansion exhibited by almost all of the populations analysed.

Post-glacial population expansion has been reported for other crustacean species within New Zealand waters (Nikula et al. 2010), and other areas of the world, including the Caribbean Sea and the North Atlantic Ocean (Maggs et al. 2008; Naro-Maciel et al. 2011). Post-glacial population expansion events among marine lobsters have largely been attributed to an increase in available habitat (Kennington et al. 2013). Lower sea-levels during the LGM do not appear to have significantly altered the amount of habitat that might have been available to NZ scampi, but may have altered habitat composition and other faunal elements. Populations of European scampi (Nephrops norvegicus) in the North Atlantic Ocean have recolonised areas after local extirpation potentially caused by glaciation (Stamatis et al. 2004; Pampoulie et al. 2011), and their genetic homogeneity has been attributed to this recent range expansion (Stamatis et al. 2004; Pampoulie et al. 2011). NZ scampi occupy similar habitat to European scampi and may have been subject to comparable extirpation events.

Management implications for New Zealand scampi

NZ scampi are currently managed as 11 separate units within the New Zealand Exclusive Economic Zone under New Zealand's quota management system. The genetically distinct populations reported in this study are currently managed separately. Photographic and trawl-based surveys of NZ scampi report stable population sizes for the areas examined within this study (Cryer et al. 2003; Tuck 2013, 2014; Tuck et al. 2015), and the estimates of genetic diversity and effective population size reported here remain high and show no indication of ongoing genetic bottleneck effects. Well managed commercial harvest of NZ scampi is unlikely to greatly impact their overall population dynamics and genetic diversity as populations are only exploited within a small proportion of their overall range. Although they live within burrows in areas of mud substrate, much of their range also incorporates rocky ground within such areas, which is not conducive to trawl-based fisheries (Wahle et al. 2011). A pot-based scampi fishery has been explored (Major et al. 2017), but does not appear to be commercially viable.

Photographic surveys of NZ scampi off the Auckland Islands observed a greater proportion of individuals emerged from their burrows compared to those of other surveyed management areas (Tuck 2015). If emergence levels for NZ scampi are greater at the Auckland Islands, their susceptibility to trawl-based fisheries is increased, as burrows offer protection from the light trawl gear employed which only catches emerged individuals (Cryer et al. 2003). NZ scampi from the Auckland Islands contain significant genetic variation unique to the Auckland Islands region, and this variation should be conserved. The current spatial arrangement of NZ scampi management areas has been designed to consider distinct NZ scampi grounds and fisheries (NZ scampi populations with similar recruitment trends and demographic parameters, and similar bycatch species) separately. This approach treats populations either side of East Cape, off the Auckland Islands, and on the Chatham Rise as separate management units, and with sustainable harvesting within these units, this continued approach should maintain stable population sizes and genetic diversity. We recommend that NZ scampi populations either side of East Cape, off the Auckland Islands, and on the Chatham Rise continue to be managed separately and are not overexploited. The present management areas and allocated harvests should achieve this.

Acknowledgements

We would like to acknowledge the following: Rick Webber (Museum of New Zealand Te Papa Tongarewa) for constructive comments and discussion regarding the life cycle of NZ scampi, and access to Te Papa's collections. Keith Probert (University of Otago) for comment on the occurrence of NZ scampi off the south-eastern South Island. Sadie Mills for access to the NIWA invertebrate collection, and the base map used in Figure 1. Nic Rawlence for constructive comments on an earlier version of this manuscript. Thank you to the crew of the RV Kaharoa and FV Bilyara for obtaining the specimens used for this study.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The DNA sequences generated for this study are available on Figshare. DOI:10.6084/m9.figshare.11465565

ORCID

Peter A. Ritchie http://orcid.org/0000-0002-8351-7931

Additional information

Funding

This work was supported by funding from Victoria University of Wellington and New Zealand’s National Institute of Water and Atmospheric Research.