Introduction

Several recent studies have focused on the issue of demographic independence, and its use in determining appropriate management units (Palsbøll et al., 2007; Waples et al., 2008; Waples and Naish 2009). The argument is that high power can be achieved with large numbers of samples and markers providing evidence for FST values significantly larger than zero, when these values may at the same time be too small to indicate demographic independence. These papers cite work by Hastings (1993) on the population dynamics of population pairs, which concludes that demographic dependence is reached when migration is >10%. It was therefore suggested that a migration level m of below 10% should be used as a criterion for assigning populations to separate management units (Waples and Gaggiotti, 2006; Palsbøll et al., 2007). However, given the large effective population sizes often estimated for marine fishes and very low associated FST values, it may not be possible to reject the null hypothesis of genetic panmixia (see below). Therefore, it is important to be able to draw on ecological and life history data to assess when apparent genetic panmixia may be consistent with demographic dependence. Long-range or distribution-wide panmixia has been reported fairly rarely (for example, Garber et al., 2005; Tolley et al., 2005; Palm et al., 2009; White et al., 2009), but was the default expectation for pelagic marine fish species until various studies demonstrated patterns of population structure, sometimes over small geographic scales (for example, Knutsen et al., 2003). In this study, we address the general question of what life history characteristics may lead to panmixia by investigating the population genetics of a deep-sea fish, the blue hake (Antimora rostrata).

Deep-sea fish species are generally characterized by slow growth rate, late maturity, low fecundity and long potential lifespan. They are increasingly threatened by commercial fishing either as the target species or as bycatch, and because of their life histories, are vulnerable to and theoretically slow to recover from overexploitation (Baker et al., 2009). The ecology of deep-sea fishes is often poorly understood. There is therefore interest in using molecular approaches to determine population structure, in order to inform management decisions, but also to make inferences about the ecology and evolution of these species.

The blue hake is a member of the gadiform family Moridae. It has a worldwide distribution, being present in all oceans from 62°S to 62°N, with the exception of the North Pacific, where it is replaced by its sister species, Antimora microlepis (Small, 1981; Turnov, 1992). A. rostrata is found at depths of 200 to 3000 m although abundance seems to peak at depths >1400 m (Fossen and Bergstad, 2006). Although the species is not commercially important, it is a common (often dominant) bycatch in deep-water fisheries (Kulka et al., 2003; Fossen and Bergstad, 2006). As a consequence, according to catch data from standardized research-trawl surveys, the blue hake is estimated to have declined in relative abundance by approximately 93% in the western North Atlantic over the period 1978–1994 (Devine et al., 2006). The blue hake from that region therefore qualifies as critically endangered according to International Union for Conservation of Nature (IUCN) criteria, as the proposed decline in abundance has been >90% over 10 years (IUCN, 2001).

Previous studies of A. rostrata have found significant differences in size between males and females, with females larger and heavier than males (Iwamoto, 1975; Wenner and Musick, 1977; Small, 1981; Gordon and Duncan, 1985; Magnússon, 2001; Kulka et al., 2003; Fossen and Bergstad, 2006). These studies also reported strong patterns of assortment by sex and by size, with larger individuals and a strong female bias in deeper sites. The nature of this distribution appears to be determined by temperature. A. rostrata occur within a narrow range of temperatures (mainly from 3 to 4 °C, but see Magnússon, 2001), and at different latitudes individuals distribute themselves differently with respect to depth, probably as a result of local slope conditions.

Magnússon (2001) reports data on maturity stages in A. rostrata from waters to the west and south of Iceland. Spawning specimens were not observed in any samples, although newly spent specimens were observed in March and July (but not in October). Maturing fish were present in July but were particularly numerous in October. Most of the maturing and newly spent specimens were observed at depths >1400 m. Magnússon therefore suggests that spawning (at least in Icelandic waters) takes place in late autumn or winter. Counts of presumed annuli from otoliths indicated a maximum age of 25 years (Fossen and Bergstad, 2006), with females appearing to live longer than males. However, because of paucity of spawning or spent individuals, it has not been possible to estimate age at maturity.

The distribution of individuals by size and sex has led several investigators to suggest that the area off Greenland might act as a nursery for the species, and that there is gradual migration to other areas (Wenner and Musick, 1977; Fossen and Bergstad, 2006). However, eggs, larvae or spawning fish have not been observed off Greenland (Kulka et al., 2003), and the occurrence of small A. rostrata along the slopes to the west of the British Isles (Gordon and Duncan, 1985) suggests that spawning might be more widely distributed. If there is a single nursery and gradual migration, we might expect to find genetic panmixia or genetic differentiation by age class. If spawning sites are more widely distributed, we might expect to find population structure between geographic samples. It is also possible that population decline and isolation by distance will have generated geographic patterns of genetic differentiation in this species. When dispersion is primarily driven by larval drift, an isolation by distance pattern could be expected, either by geographic distance or along current paths (for example, Knutsen et al., 2007). In this paper, we used 13 microsatellite loci to test the hypothesis that the proposed single dominant spawning site and subsequent gradual migration, or some level of periodic re-distribution during reproduction (necessary if adult fish remain geographically segregated by sex at other times) has led to genetic panmixia across an extended geographic range. Given an estimated decline of 93% for blue hake in the western North Atlantic, we further consider a simple model with general relevance, and the implications for the rate of migration that could lead to demographic dependence.

Materials and methods

Sampling

In total, 393 individuals of A. rostrata were obtained from 11 sample sites in the North Atlantic (Table 1 and Figure 1). These samples were obtained from the mid-Atlantic ridge (MAR) and from North of the Azores Archipelago from the MARECO expedition in July 2004 and two ECOMAR expeditions in July 2007 and 2009. From the Porcupine Seabight, samples were obtained from the RRS Discovery Cruise no. 252 in April 2001. From West Greenland, samples were obtained from Ole Jørgensen from a research cruise in October 2008.

Table 1 Sampling information and summary statistics for each sample site
Full size table
Figure 1
figure 1

Sample sites in the North Atlantic. Gray scale shading denotes height above sea level in meters. Dashed line shows approximate location of the sub-polar front. The top left of the map is 70°N, 70°W, and the bottom right is 40°N, 10°E.

Full size image

Distribution by sex and length

For all sites except Porcupine we had biometric data, and were able to test for geographic distribution by sex and length. All analyses were carried out using the R statistical language (R Development Core Team, 2009).

DNA extraction and genotyping

DNA was extracted from specimens using a phenol–chloroform protocol (after Hoelzel, 1998). A total of 17 microsatellite loci were PCR amplified. These were previously described by White et al. (in press). PCR reactions were carried out using Qiagen multiplex PCR kits (Qiagen, Crawley, UK) in a total volume of 10 μl per reaction. This volume comprised 5 μl of multiplex PCR master mix, 4 μl of primer master mix (with each primer at a concentration of 0.5 μM) and 1 μl of template DNA. PCR was performed with the profile: hot-start activation at 95 °C for 15 min; 35 cycles of 30 s denaturing at 94 °C; 90 s annealing at 57 °C; 60 s extension at 72 °C; final extension at 60 °C for 10 min. Loci were multiplexed in three groups: group 1 (ARos1, ARos2, ARos4, ARos6, ARos8), group 2 (ARos12, ARos14, ARos20, ARos21, ARos22) and group 3 (ARos29, ARos32, ARos33, ARos37, ARos38, ARos39, ARos40). Forward primers were 5′ fluorescently labelled using either FAM, HEX or TAMRA (Eurofin, Ebersberg, Germany) for screening of alleles on the ABI 3730 sequencer. Alleles were scored using PEAK SCANNER v1.0 (Applied Biosystems, Warrington, UK) with GS 500 as the internal size standard. Alleles were binned using the program FLEXIBIN v2 (Amos et al., 2007).

Genetic analysis by sample site

Linkage disequilibrium between loci was tested using GENEPOP 3.4 (Raymond and Rousset, 1995). For each sample site, expected heterozygosities, allelic richness and FIS were calculated using GENEPOP. Significance of deviation from Hardy–Weinberg equilibrium was evaluated using exact tests. The program MICROCHECKER v2.2.3 (Van Oosterhout et al., 2004) was used to test for allele scoring errors and the existence of null alleles. MICROCHECKER found evidence of null alleles in more than one sample for the loci ARos8, ARos22, ARos32 and ARos37. When these loci were included in our analyses, FST between some pairs of populations increased, but this occurred in populations of small sample sizes and was clearly influenced by the presence of null alleles (data not shown). Therefore, further analyses were conducted without these four loci.

Including loci under selection can sometimes distort population genetic analyses (White et al., 2010). Therefore, we used a Bayesian approach to detect loci under selection, implemented in the program BAYESCAN (Foll and Gaggiotti, 2008). FST values reflect contributions from locus-specific effects, such as selection, and population-specific effects such as genetic drift and immigration rates (Balding and Nichols, 1995). BAYESCAN uses a hierarchical Bayesian approach to estimate the locus and population effects on these FST values. Following the suggestion of Foll and Gaggiotti (2008), we used the ‘decisive’ criterion under Jeffrey's scale of evidence to identify outlier loci, in order to minimize the false-positive rate and to maximize the true-positive rate. No loci emerged as outliers under selection.

We used the program POWSIM 4.1 (Ryman and Palm, 2006) to assess the statistical power of our analysis to reject the null hypothesis of genetic homogeneity for different combinations of number of samples, sample sizes, number of loci, number of alleles and allele frequencies for a hypothetical degree of true differentiation (quantified as FST).

Overall and pairwise FST (Weir and Cockerham, 1984) among sample sites in the North Atlantic was calculated using ARLEQUIN 3.1 (Excoffier et al., 2005), and the significance of the difference between the observed FST and zero was evaluated using exact tests. Fisher's exact test (Ryman et al., 2006) in GENEPOP was used to test for significant differences in allele frequencies between populations. As some investigators have questioned the use of FST as a measure of genetic differentiation (Hedrick, 2005; Jost, 2008), we also calculated Jost's D over all sample sites and for pairwise comparisons using SMOGD (Crawford, 2010).

FST values were transformed into FST/(1–FST) and regressed onto the log distance between population pairs. Significance was assessed using a Mantel test with 10 000 permutations, implemented in GENEPOP.

We tested the effect of different environmental factors on the genetic structure of roundnose grenadier using a hierarchical Bayesian method, implemented in the program GESTE (Foll and Gaggiotti, 2006). This method estimates population-specific FST values and then relates them to environmental factors using a generalized linear model. The environmental variables considered here were depth, latitude and longitude, year of sampling, side of ridge (west vs east) and side of sub-polar front (north vs south) generating 64 different models with all combinations of factors (including a constant term). GESTE was run using default parameter values for the Marcov Chain Monte Carlo (MCMC) analysis. Posterior probabilities of including a given factor in a model were obtained by summing the posterior probabilities of all the models that include this factor.

We also conducted analysis of molecular variance in ARLEQUIN 3.1 (Excoffier et al., 2005), with several different hierarchical structures relating to barriers that are thought or known to influence gene flow (White et al., 2010). First, sample sites were grouped according to the side of the sub-polar front on which they fell (north vs south). Second, they were grouped according to the side of the ridge on (west vs east). Third, sample sites were placed in one of four groups, based on side of the sub-polar front and side of ridge. Finally, sample sites were grouped into six broad geographic clusters (see Table 1).

Genetic analysis by age class

Total length data were available for all samples except for those from the Porcupine region. Individuals were assigned to approximate age classes on the basis of their total length, using the Von Bertalanffy Growth Function parameters given in Fossen and Bergstad (2006) with the formula Age=ln((L∞–L(t)/(L∞–L0))/–k, where k is the growth parameter, L0 is the length at age 0, L(t) is the length at age t and L∞ is the ultimate length attained at age infinity. Although females ultimately grow to be larger than males, the Von Bertalanffy curves of males and females do not differ significantly (Fossen and Bergstad, 2006). Therefore, we used the parameters estimated using combined data from males and females. Age classes were then adjusted according to the year of sampling. We then followed the approach of Knutsen et al. (2009) in comparing the effect of age classes on genetic differentiation. Genetic differences aij (Rousset, 2000) were calculated between pairs of individuals (ij), averaged over the 13 loci. The measure aij was then used as the dependent variable in a generalized additive model:

where tij is the estimated difference in age between two individuals, and cij is a Boolean variable indicating whether the two individuals were caught at the same location or not. In the above formula s(x) are smoothing spline functions with degrees of freedom chosen by cross-validation, and allowing for the expected nonlinear effects of age differences when generations overlap (Jorde and Ryman, 1996). The model was fitted using the generalized additive model library in the R statistical package. Significance testing was carried out by bootstrapping of individual covariates to generate null distributions.

Individuals were also clustered on the basis of age to perform population-based analyses. Age classes with <10 individuals were excluded. Again, overall and pairwise FST (Weir and Cockerham, 1984) values among age classes were calculated using ARLEQUIN, and the significance of the difference between observed FST and zero was evaluated using exact tests. Fisher's exact test in GENEPOP was used to test for significant differences in allele frequencies between populations. Pairwise FST values then regressed onto the temporal distance between age classes. Significance was assessed using a Mantel test with 10 000 permutations, implemented in GENEPOP.

Genetic analysis by sex

Individuals from all localities except Porcupine were pooled into two groups, males and females. FST between males and females was calculated using GENEPOP, and the significance of the difference between observed FST and zero was evaluated using an exact test in ARLEQUIN and using Fisher's exact test in GENEPOP.

FST, dispersal, demographic independence and power

We can ask how much power would be required in order to detect demographic independence. Hastings (1993) suggests that populations become demographically dependent when m=0.1. Takahata (1983) shows that in a finite island model with d demes

Assuming that two distant sample sites could be considered separate demes exchanging migrants, we estimated the FST value that would be generated given demographic independence (m<0.1). We chose Mar6 and Greenland for this analysis. These sites were separated by 2095 km and had large sample sizes. Ne values (see below) for these sample sites were used with Takahata's formula to estimate FST, and we assessed our power to detect such an FST using POWSIM.

Given data suggesting ∼90% depletion of blue hake stocks in the western North Atlantic over a period of ∼20 years, we devised a simulation model that tested the migration rate required to equalize the size of two populations (A and B), one depleted and one not, given different parameters associated with intrinsic growth (r) and environmental stochasticity (z). Depletion was based on a harvest rate required to reduce the population from a hypothetical census size of 250 000 to 25 000 over 20 years, given a specific intrinsic rate of growth. Density-dependent growth is given by:

where r is the intrinsic population growth rate and Kt is the carrying capacity. Harvesting (hA or hB) occurred at some set rate for each population and modified the population size so that A′′=A′(1–hA). Finally, we assumed that migration was represented by some set proportion, m, of each population moving to the other population each year. Thus, the following year's population size for population A was given by:

We are particularly interested in the effect of migration when only one of the populations is harvested. For this reason, we set hA=0. Estimates of r for similar marine fish vary but are unlikely greatly to exceed r=1 (Myers et al., 1997; Baker et al., 2009). To determine the effect of r on the outcomes of our simulations, we used values of r of 0.1, 0.5 and 0.9.

Simulations were affected by environmental stochasticity. Stochasticity was incorporated by varying Kt, the carrying capacity, from year to year. Specifically, Kt=ztK̄, where zt is a normally distributed random variable with a mean of 1.0 and a s.d. of σ, and K̄ is the mean carrying capacity. Carrying capacities of the two populations were assumed to vary independently, representing different conditions in different parts of the species' range. We assessed the impacts of a variety of different levels of stochasticity (σ=0–0.6) but, because qualitative outcomes were unaffected by the level of stochasticity, we present only the case with a conservative level of stochasticity (σ=0.1).

We assessed the outcomes of simulations in two ways. First, we looked at the effect on both populations of starting with 250 000 fish in each, and then harvesting one population only for 20 years (that is, hA=0, hB >0). The harvesting rate was that which would lead to a 90% deterministic decline of an isolated fish stock harvested at that rate (that is, with constant effort) for 20 years (without compensating for the increasing migration rate). This result was expressed by plotting, as a function of m, the mean size of both populations at the end of a 20-year period. Results were averaged over 5000 replicate simulations. Our second approach was to simulate the sizes of the two populations under the same harvesting regimes but over a longer period. Specifically, we ran simulations for 110 years and assessed the extent to which annual sizes of the two populations were correlated over the last 100 years of the simulation. To determine the mean level of correlation associated with any parameter set, we repeated this process 1000 times. All simulations were conducted in R.

Effective population size and population size changes

The software DIYABC (Cornuet et al., 2008) was used to estimate the effective population size for Mar6 and Greenland. We also estimated effective population size assuming that there is only one population in the whole of the North Atlantic. For all three runs of DIYABC, we defined two models, one assuming constant population size and one allowing a population size change at t generations in the past. Previous uniform distributions were set for Ne (lower limit 100 and upper limit 100 000), t (lower limit 1, upper limit 10 000), mean μ (lower limit 1 × 10−5 and upper limit 1 × 10−3) and a gamma distribution for locus μ (lower limit 1 × 10−5 and upper limit 1 × 10−2, with shape 2.0). The stepwise mutation model was used and one million replicate runs were performed to generate the reference tables.

We also used Garza and Williamson's M ratio to detect changes in population size. M ratios were calculated using a Perl program. Expected values of M under demographic equilibrium were calculated using the program Critical_M (Garza and Williamson, 2001). Following the recommendations of the authors, 12% of mutations were allowed to be greater than one step, with a mean size of larger mutations of 2.8 steps. Theta (4Neμ) values were approximated using the DIYABC output.

The software Bottleneck (Cornuet and Luikart, 1997) was used to test for recent reductions in population size. We used both the single-step model and the two-phase model (88% of single-step mutations and a variance of 30) of microsatellite evolution, and used the sign test and the Wilcoxon test to determine whether any of the sample sites showed an excess of heterozygosity, as expected after a severe bottleneck.

Results

Distribution by sex and length

As in previous studies (which used overlapping data sets), we found a strong female bias in all sites except Greenland. From the Greenland site, there were 34 females (37%) and 58 males (63%), whereas over all other sites there were only 5 males (1.6%—no biometric data were available for the Porcupine site). The sexes differed markedly in total length (the female mean was 495 mm, while the male mean was 341 mm: t=−11.31, P<0.001). However, most of this effect may be due to location of sampling, as in Greenland both small males and females were found. When location and sex were included in an analysis of variance, only location (F=83.82, d.f.=9, P<0.001) and the interaction term (F=5.18, d.f.=4, P<0.001) emerged as significant, while the main effect of sex was not significant (F=1.49, d.f.=1, P=0.223). The effect of location appears to be mediated by latitude and depth of sampling, as individuals from Greenland were small, individuals from the North of the Azores Archipelago tended to be large, and individuals from elsewhere on the MAR were of intermediate length.

Genetic analysis by sample site

None of the loci used in this study showed any evidence of linkage disequilibrium. Basic diversity statistics for each sample site are shown in Table 1, and genetic diversity was high in all samples. FIS values were mostly positive, and four remained significant after Bonferroni correction (Table 1). FST over all populations was not significantly different from zero (FST=0.001; ARLEQUIN exact test: P=0.800; Fisher's exact test: P=0.11). None of the pairwise FST values were significantly different from zero following exact tests in ARLEQUIN. Three pairwise comparisons showed significant differences in allele frequencies at the 5% level using Fisher's exact test, but not after Bonferroni correction was applied (Table 2). We also calculated q-values as a less conservative alternative to account for type 1 errors (Storey 2002), and no significant values were found (data not shown). Pairwise D values ranged from −0.011 to 0.016. Overall D was 0.002 (a low value comparable to that seen for FST).

Table 2 Genetic differentiation between sample sites
Full size table

POWSIM suggested that with our loci and our sampling regime, we would have 99.5% power to detect significant genetic differentiation when FST=0.002, 90% power when FST=0.0015, 66.5% power when FST=0.001 and 30% power when FST=0.0005.

There was no relationship between FST /(1–FST) and the log distance between sample sites (r=0.266, P=0.20. Mantel test, 10 000 permutations). Furthermore, none of the analysis of molecular variances revealed significant genetic differentiation between groups (see Materials and methods). In each case, the vast majority of the variance was between individuals within populations (Table 3). From the GESTE analysis, no environmental factor emerged as a significant predictor of genetic differentiation between samples (results not shown). In each case, only the model including only the constant term received any support.

Table 3 Results of AMOVAs
Full size table

Genetic analysis by age class and sex

The results of the generalized additive model analysis are shown in Table 4. There was a significant relationship between temporal distance and genetic distance, although the effect size was extremely small. Individuals caught at the same location also appeared to be significantly more genetically differentiated than those caught at different locations (positive effect of the covariate c). When individuals were clustered into age classes, FST over all age classes was 0.000 and not significant. There was no relationship between the pairwise genetic distance and temporal distance between age classes (r=−0.131, P=0.77. Mantel test, 10 000 permutations). POWSIM suggested that with our loci and the numbers in each age class, we would have 90% power to detect significant genetic differentiation between age classes when FST=0.002, 75% power when FST=0.0015, 48.5% power when FST=0.001 and 26.5% power when FST=0.0005. FST between males and females was not significantly different from zero (FST=–0.0007; exact test: P=0.14).

Table 4 Results of GAM model testing for the effects of sampling location and age class on genetic differentiation among A. rostrata individuals
Full size table

Effective population size and population size changes

DIYABC was used to estimate effective population size and changes in population size over time. For all three data sets, direct estimates and logistic regression estimates gave much greater support to the model allowing population size to vary. For Mar6, the posterior median of current Ne was estimated to be 14 500 (95% confidence limits: 5660, 36 600), for a median mutation rate of 6.71 × 10−4 (2.62 × 10−4, 9.85 × 10−4). The posterior distributions for past Ne and t, the time of population size change, were broad, with past Ne having a median of 75 600 (95% confidence limits: 35 100, 98 900), and t having a median of 5600 generations (95% confidence limits: 839, 9 780). For Greenland, the posterior median of current Ne was estimated to be 14 500 (95% confidence limits: 6020, 35 900), for a median mutation rate of 7.16 × 10−4 (2.96 × 10−4, 9.87 × 10−4). The posterior distributions for past Ne had a median of 81 500 (95% confidence limits: 42 800, 99 200), and t had a median of 4880 generations (95% confidence limits: 683, 9710).

Considering all the sample sites as one population, the posterior median of current Ne was estimated to be 15 700 (95% confidence limits: 5960, 43 500), for a median mutation rate of 6.96 × 10−4 (2.71 × 10−4, 9.89 × 10−4). The posterior distributions for past Ne had a median of 77 600 (95% confidence limits: 33 500, 98 900), and t had a median of 6740 generations (95% confidence limits: 547, 5690).

These estimates of Ne and mutation rate were used to calculate theta (4Neμ), to determine the expected value of M under demographic equilibrium using the program Critical_M. With theta set to 40, approximating DIYABC's estimate of the current population size, the critical M value was 0.77 for Greenland and 0.74 for Mar6. With theta set to 210, representing the historical population size, the critical value was 0.70 for Greenland and 0.65 for Mar6. Although M values were less than these critical values for the lower estimates of Ne (Table 1), they were higher than those estimates using the historical Ne. Given the uncertainties involved, there is no clear evidence of a recent genetic bottleneck in any population.

None of the sample sites showed evidence of heterozygosity excess under either the two-phase model or the stepwise model of microsatellite evolution (results not shown), providing no genetic evidence for a recent bottleneck.

FST, dispersal, demographic independence and power

Using Takahata's formula for the two demes Mar6 and Greenland, with our estimates of current Ne set at 15 000, and with d=2, we get an estimate for FST≈0.000042 when m=0.1 (Hastings, 1993). This value is much lower than could be detected using our data set. This level of FST was also lower than the minimum value that could be assessed using POWSIM. However, running POWSIM at the limits with 50 loci (sampled randomly several times from the 13 loci used here), we obtained a power of 7% to detect FST=0.0001.

Outcomes from our simulations showed the expected (trivial) result for the two extremes (when m=0 or 0.5), but a non-linear relationship between population size and intermediate migration rate (Figure 2). Furthermore, the point at which the two populations converge in size was dependent on the intrinsic rate of growth, and in all cases greater than m=0.1 (Figures 2a–c). The extent to which the two populations are correlated (with respect to annual changes in population size, rather than overall size) also varied with their intrinsic growth rate, as well as with the extent of migration (Figures 2d–f). In this case 100% correlation was achieved by m=0.1 for the lowest value of r, but not until m=0.5 for the higher values of r.

Figure 2
figure 2

The consequences of migration for the fates of two populations, one of which is harvested. (a–c) Final population sizes following a 20-year simulation in which one population (broken line) is harvested and the other (solid line) is not (lines show average over 5000 Monte Carlo simulations). (d–f) Mean correlation (Pearson's r) between the two populations over 100 years of simulation (averaged over 1000 Monte Carlo simulations). See Materials and methods for further details. K̄=250 000, σ=0.1 in all simulations. Intrinsic rate of population growth, r, and rate at which population (b) was harvested, hB, increase from left to right: (a, d) r=0.1, hB=0.16; (b, e) r=0.5, hB=0.38; (c, f) r=0.9, hB=0.56. Shaded areas show 95% confidence intervals for simulation outcomes.

Full size image

Discussion

We found no evidence of spatial genetic structure in A. rostrata, despite having sufficient power to reject panmixia when FST is greater than approximately 0.0015. As our FST could not be distinguished from zero, we accept the hypothesis of panmixia and therefore demographic dependence across the sample range. However, we also used a simple simulation model to consider further the conditions under which demographic dependence becomes likely. Our model illustrates that demographic correspondence is highly dependent on the intrinsic growth rate of the population. Baker et al. (2009) give estimates of r for two deep-sea fishes, Macrourus berglax and Coryphaenoides rupestris, ranging between 0.082 and 0.568, while Myers et al. (1997) gave estimates of r between 0.15 and 1.03 for the Atlantic cod, Gadus morhua. On the basis of our simulations, convergence of population size may occur at relatively low rates of migration when r is low, but at much higher rates when r is high. For a credible range of r for a marine fish (between 0.5 and 0.9) it seems that demographic correspondence would occur after the migration rate was considerably higher than 10% (closer to 40%, see Figure 2, although the depleted population starts to impact the non-harvested population even at low migration rates). At the same time, we show that it would require many more loci (perhaps 100 s of loci) to provide the resolution to determine an FST value that reflects separate populations sustaining such high levels of migration, and such low FST values are unlikely to be biologically meaningful. On the basis of populations with Ne comparable to that estimated here and with loci showing levels of diversity comparable to those reported here, we show that many more than 50 loci would be required even to assess an FST resulting from a migration rate of 10%. There is, therefore, a considerable gap between the migration rate level at which we have the potential to reject panmixia, and the migration rate that may be necessary to generate demographic correspondence in this species.

Our estimate of current Ne for A. rostrata for the North Atlantic as a whole was high, at approximately 15 000. Despite recent studies suggesting a massive population decline in this species (Devine et al., 2006), and our analyses suggesting a long-term decline in Ne, our results suggest that the population as a whole may still be large. Recently, Cushman et al. (2009) found genetic homogeneity between different populations of the gag (Mycteroperca microlepis) from the Gulf of Mexico and northwest Atlantic. In that study, Ne estimates were also very high (≈16 500). The authors attempted to use a coalescent-based simulation approach to infer the level of migration that would lead to such low levels of genetic structure; however, they were unable to do so accurately. When FST is so low, these approaches are expected to fail, because of the non-linear relationship between Nm and FST. Within the margins of error, a low value of FST is consistent with Nm ranging from the tens to infinity (Whitlock and McCauley, 1999).

Other recent studies using microsatellite loci have found populations in the North Atlantic to be genetically panmictic (Palm et al., 2009; White et al., 2009). For marine species occupying a large range of environments and with large effective population sizes, neutral loci may exhibit little divergence despite fairly limited gene flow, whereas loci under strong selection may display a high degree of differentiation (White et al., 2010). Hence, using loci under selection could aid in identifying ecologically and genetically relevant units, while the distribution of different alleles in time and space could also be used to infer dispersal parameters (Barton and Gale, 1993; Hare et al., 2005).

However, using the genetic data at our disposal, we cannot determine whether the blue hake in the North Atlantic is really one demographic unit. Nevertheless, the pattern of distribution by size and sex is consistent with this interpretation. We found a non-random distribution of individuals in relation to size and sex, and relatively high FIS levels within most sample-sets. This suggests that standard population comparisons may not be appropriate in this species and that some kind of temporal Wahlund effect may be operating (Hoarau et al., 2005), for example, because of sweepstake effects (in which most surviving larvae are from a relatively small number of females) every breeding season (Hedgecock, 1994). In common with previous studies of A. rostrata, we found that samples from off Greenland contained both males and females, whereas samples from elsewhere consisted almost entirely of females. Individuals from Greenland also tended to be smaller than those from elsewhere. Breeding individuals of A. rostrata have rarely been reported (off Iceland, Magnússon, 1998), and eggs and larvae have never been found (Kulka et al., 2003). This may reflect the difficulty of identifying juvenile fish to the correct species (Peter Rask Møller, personal communication). As the smallest specimens of A. rostrata so far described have been found off Greenland, it is tempting to regard this region as a breeding ground and/or nursery. Nevertheless, we are far from having a clear picture of the distribution of this species and only in a few places has the full-depth range been sampled (Fossen and Bergstad, 2006). Males and small females may occur in large numbers on the MAR at depths shallower than those sampled. Similarly, large females may occur in greater numbers off Greenland at depths greater than those sampled.

Given the sex and size distribution of this species, and lack of genetic differentiation between localities, we assigned individuals to age classes on the basis of length data and year of capture. Taking age classes to represent populations revealed no significant structure between years and the positive FIS values remained. However, further analyses revealed statistically significant but very small effects of age class and area of sampling on the genetic differentiation between individuals. Individuals from the same geographic locality were marginally (but significantly) more genetically differentiated than those from different localities. This is opposite to the usual prediction, but may be indicative of an extremely high exchange of individuals between localities, bringing together individuals that are less related than expected by chance. Such mixing is consistent with the spatial structuring of the population with regard to size and sex, and the lack of genetic differentiation between males and females, despite their rather extreme geographic segregation. As a working hypothesis, we suggest that the population of A. rostrata within the North Atlantic is one panmictic unit, consistent with the genetic data. At some point males must come into contact with mature females, and the sequential re-distribution of sexes would lead to mixing. Different age classes may exhibit a small degree of genetic differentiation because of drift. However, this differentiation may then be lost by mixing of individuals on an ocean-basin scale, followed by segregation by depth according to age and sex. Only by identifying breeding areas and conducting genetic analysis on juveniles would it be possible to test this hypothesis fully (Burford and Larson, 2007). Other marine species with similar life history may also show large-scale panmixia because of the periodic redistribution of the sexes for mating.