Levels of neutral genetic diversity largely reflect effective population size (Ne), which is influenced by physical and biological environmental variables. While large populations of marine fishes generally harbor higher diversities than freshwater species, historical demography or the effects of natural selection may significantly reduce Ne and differentially affect genomic diversities. Here, we surveyed levels of genetic diversity and examined genetic structure among populations of the Atka mackerel Pleurogrammus monopterygius across most of its geographic range by screening variation at nine nuclear microsatellite DNA markers (n = 745) and in a 468-base-pair segment of the mitochondrial DNA (mtDNA) control region (n = 119). Samples from Japan to the western Gulf of Alaska were collected between 2004 and 2006 at six locations, including temporal replicates at two sites. Microsatellite allele frequency homogeneity across the North Pacific indicated an apparent lack of genetic population substructure. While levels of polymorphism at microsatellite loci were typical for marine fishes (haplotype diversity h = 0.34–0.96), mtDNA control region diversity was extremely low (nucleotide diversity = 0.00007; h = 0.033). Only three mtDNA haplotypes, two occurring as singletons, were detected among 119 individuals. The strong contrast between microsatellite and mtDNA diversities appears to be due to the smaller Ne for mtDNA, perhaps resulting from population bottlenecks during postglacial colonizations of the central North Pacific or the effects of natural selection on mtDNA.
Gene diversities of neutral molecular markers are most influenced by effective population size (Ne; Crow and Kimura 1970). Diversities are expected to be high in large, stable populations because the magnitude of random drift is less, leading to the retention of a larger proportion of existing diversity and a greater number of new mutations. Other factors also influence Ne and, hence, gene diversities. One important factor is temporal variability in population size; populations that experience bottlenecks (Nei et al. 1975) or episodic extinctions and recolonizations (metapopulation effects; Hanski and Gilpin 1991) can have low levels of diversity despite their large contemporary population sizes. A skewed sex ratio of reproducing individuals can also reduce Ne (Crow and Kimura 1970), and the mode of inheritance can influence the gene diversity of a molecular marker. Biparentally inherited nuclear genes are expected to have Ne values that are four times those of uniparentally inherited genes encoded by mitochondrial DNA (mtDNA; Birky et al. 1989); nuclear genes are thus less susceptible to erosion of neutral genetic variation via random drift during founder events or population bottlenecks. Hence, at some values of Ne, random drift can reduce mtDNA gene diversities without affecting nuclear gene diversities (Grant and Leslie 1993).
In the marine realm, high levels of adult mixing or larval drift (gene flow) tie populations together such that marine species are expected to have large Ne values (Waples 1998). This predicted broad-scale panmixia has been reported for many open-ocean species (Palumbi and Wilson 1990; Dannewitz et al. 2005; but see Wirth and Bernatchez 2001; Gilbert-Horvath et al. 2006). Another prediction is that high gene diversities should reflect large contemporary population sizes in marine species. In general, marine fishes tend to harbor greater allozyme and microsatellite diversities than freshwater fishes, which are found in smaller, spatially partitioned habitats that limit Ne (Ward et al. 1994; DeWoody and Avise 2000). Gene diversities are maximal at migration–drift equilibrium, and the rate of approach to this maximum is inversely related to Ne (Nei et al. 1975). However, the expected high levels of diversity in abundant marine species may be altered by periodic disturbances—generally associated with climatic cycles—that prevent populations from achieving this equilibrium.
In this study, we estimated genetic diversity and examined the genetic population structure of North Pacific Atka mackerel Pleurogrammus monopterygius (Family Hexagrammidae). This species is distributed along continental shelves and slopes to 500 m from Japan, along the eastern Kamchatka Peninsula, through the Aleutian Islands archipelago (the region of highest abundance), to the Gulf of Alaska and southeast Alaska (Allen and Smith 1988; Lauth et al. 2007ba). Females are obligate demersal spawners that deposit egg masses on rocky substrates in areas of moderate to strong currents, and the egg masses are cared for and defended by the males (Lauth et al. 2007ba). Numerical sex ratios are approximately equal (S.A.L., unpublished data), and Atka mackerel exhibit a polygynandrous mating system (Canino et al. 2010). Nesting begins in June, and the overlapping mating and brooding seasons may last up to 7 months (Lauth et al. 2007b). Newly hatched larvae are neustonic (Kendall and Dunn 1985) and develop in pelagic habitats for 1 year before recruiting into demersal populations as juveniles. Tagging results indicate that adults form local aggregations and generally move less than 70 km (McDermott et al. 2005).
Contrasting patterns of geographic variability are evident in previous studies of growth and life history characteristics and genetic molecular markers. Phenotypic differences among areas indicate the possibility of discrete Atka mackerel populations (Lee 1985; Lowe et al. 1998). In particular, length and weight at age differ significantly among locations along the Aleutian Islands (Kimura and Ronholt 1988), and length at age varies clinally, increasing from west to east along the Aleutian Islands archipelago (Lowe et al. 1998). These differences are consistent with limited adult migration but may not necessarily reflect genetic differences that have accumulated among populations. Indeed, a pattern of geographical heterogeneity did not appear in the frequencies of allozymes (Lowe et al. 1998). A survey of 22 polymorphic allozyme loci in four samples (n = 329) along the Aleutian Islands archipelago failed to show significant frequency differences (genetic differentiation index FST = 0.004). Additionally, the level of genetic diversity (average heterozygosity H = 0.137) was exceedingly large, more than double the average reported by Ward et al. (1994) for 57 marine fish species.
The previous genetic results suggest that high levels of gene flow between populations may result from larval drift in currents, if not from the movements of sedentary adults. In this case, phenotypic differences among adults in different areas may reflect short-term responses to habitat differences in temperature and food availability. Alternatively, the observed phenotypic differences may have a genetic basis uncoupled from allozyme-encoding genes, or allozyme markers may not provide sufficient population resolution. The use of molecular markers with higher mutation rates, such as mtDNA and microsatellite DNA, may provide the resolution needed to detect population differences on the same geographic scale as phenotypic differences. The high level of allozyme diversity reported by Lowe et al. (1998) further suggests that values of Ne for Atka mackerel are large and that populations have reached mutation–drift equilibrium. However, a long history of large population sizes in the North Pacific is unlikely due to the strong ocean climate oscillations that occurred during the Pleistocene (Mann and Peteet 1994; Mann and Hamilton 1995).
The primary goal of this study was to use microsatellite and mtDNA markers to extend the genetic analysis of Atka mackerel populations. We first tested a null hypothesis of no regional differentiation among samples, and we used microsatellite DNA and mtDNA to contrast with previous results for allozymes. In comparison with allozymes, these two marker classes have been shown to provide greater resolution of the weak structure typical of marine fish species (Bentzen et al. 1996; McPherson et al. 2001; Withler et al. 2001; Cunningham et al. 2009). While both marker types showed broad-scale regional homogeneity, an unexpected result was a strong discordance between nuclear and mitogenomic diversities. Atka mackerel populations appear to be recovering from the effects of population bottlenecks, postglacial founder events, or selective sweeps on mtDNA that were sufficiently large to greatly erode mtDNA variability but not large enough to significantly reduce the diversities of nuclear genes.
Methods
Sample collection, DNA extraction, and polymerase chain reaction and sequencing protocols
Fin clips were collected from adults captured by trawling at each of four locations along the Aleutian Islands archipelago during the 2004 summer spawning season. Two of these locations were re-sampled in 2006 (Table 1; Figure 1). Additional samples originated from the Sea of Japan in 2004 and from the western Gulf of Alaska in 2005. All samples were stored in 95% nondenatured ethanol at room temperature prior to DNA extraction.
Genomic DNA was extracted using Qiagen DNeasy 96 tissue kits according to the manufacturer's instructions (Qiagen, Inc., Valencia, California). Polymerase chain reaction (PCR) was used to amplify seven microsatellite loci isolated from Atka mackerel (Pmo69, Pmo70, Pmo152, Pmo164, Pmo268, Pmo367, and Pmo399; Spies et al. 2005) and two additional loci (Oel42 and Oel32; S. Young, Washington Department of Fish and Wildlife, personal communication) isolated from a related hexagrammid species, the lingcod Ophiodon elongatus. Genotypes were determined with a 4200 LI-COR DNA analysis system (LI-COR Biotechnology, Lincoln, Nebraska) and LI-COR SAGAGT genotyping software using a double-blind scoring protocol to minimize genotyping errors.
A 468-base-pair (bp) segment of mtDNA control region (CR) was sequenced in 119 individuals from three locations spanning the sampling range (Table 1). The mtDNA CR was amplified by PCR with primers ProL and H16498 via protocols given by Withler et al. (2004). The PCR amplicons were diluted with distilled water to 5.5 ng/μL and were purified using ExoSAP-IT (USB Corp., Cleveland, Ohio). Forward and reverse DNA strands were sequenced with the same primers following standard procedures at the High-Throughput Genomics Unit (University of Washington, Seattle). Sequence contigs were assembled using Sequencher version 4.8 (Gene Codes Corp., Ann Arbor, Michigan), and aligned consensus sequences were created using BioEdit version 7.0.9.0 (Hall 1999). Sequences were deposited in GenBank (accession numbers FJ858207–FJ858209).
Microsatellite analysis
Micro-Checker version 2.2.3 (van Oosterhout et al. 2004) was used to check for null alleles, scoring errors, stuttering, and large-allele dropout. The genotypic frequency fit to Hardy–Weinberg (HWE) expectations and linkage disequilibrium were evaluated by exact tests implemented in GENEPOP version 1.2 (Raymond and Rousset 1995) using 500 Markov-chain batches of 5,000 steps/batch after 5,000 dememorization steps. Where appropriate, significance levels for multiple tests were corrected to produce a type I error level (α) of 0.05 with the sequential Bonferroni correction (Rice 1989). Global tests of heterozygote excesses and deficits by locus and across loci were made with GENEPOP. Observed H (HO), expected H (HE), and allelic richness were estimated with FSTAT version 2.9.3.2 (Goudet 2001). Heterogeneity among samples was determined using exact tests in GENEPOP with Markov-chain Monte Carlo (MCMC) parameters specified above and log-likelihood g-tests implemented in FSTAT with 10,000 data randomizations. Weir's (1996) estimator of FST for each locus and across loci was calculated between sample pairs and for all samples using FreeNA (Chapuis and Estoup 2007), and 95% confidence intervals were obtained with 10,000 bootstrap resamplings. A standardized measure of FST (G′ST; Hedrick 2005) that accounts for the effects of locus polymorphism on estimates of differentiation was calculated for each locus and over all loci.
Several approaches were used to assess geographic variation. Evidence for isolation by distance was examined with Mantel's (1967) test, as implemented in the ISOLDE routine in GENEPOP, using a matrix of FST values and geographic distances calculated as nearest distances between sampling sites with Google Earth ( earth.google.com/). Analysis of molecular variation (AMOVA) implemented in ARLEQUIN version 3.11 (Excoffier et al. 2005) was used to partition genetic variation across the geographic sample range into its spatial and temporal components. Finally, the program STRUCTURE version 2.2 (Pritchard et al. 2000) was used to infer whether samples formed geographically distinct clusters within the sampled range. The mixed-ancestry model with correlated allele frequencies was used to infer the number of populations or clusters (K) represented in the samples. Posterior probabilities for specified K-values of 1–6 were calculated with 100,000 MCMC iterations after burn-in discard of the first 50,000 iterations.
Historical signals of demographic shifts were assessed in two ways. First, allele frequency distributions of were examined with BOTTLENECK (Cornuet and Luikart 1996) to test for a recent reduction in Ne. Bottleneck events produce deficits of low-frequency alleles through drift (heterozygosity excess, distinguished from Hardy–Weinberg heterozygote excess) in comparison with the level of genetic diversity expected in constant-size populations at equilibrium (Luikart and Cornuet 1998). A single-stepwise mutation (SSM) model and a two-phase mutation (TPM) model (Di Rienzo et al. 1994) of microsatellite evolution were used, with 5% and 10% multistep mutations for the TPM simulations; significances were determined with Wilcoxon's signed rank tests (Luikart and Cornuet 1998). Second, the shapes of the microsatellite DNA trees were examined with k- and g-tests (Reich et al. 1999) as implemented by the spreadsheet macro of Bilgin (2007). The k-statistic examines individual locus trees to distinguish between (1) a deeply bifurcated tree with clusters of microsatellite alleles of similar lengths, as expected in populations of constant size, and (2) a shallower tree with nodes of most branches dating to the time of population expansion. The g-test examines the distributions of coalescences among multilocus trees. With constant population size, coalescences at the base of the trees vary considerably; however, in an expanding population, coalescences tend to have similar dates reflecting the expansion event. Significance of the k-test values was determined by the proportion of loci giving a negative k-value using a one-tailed binomial distribution, and the probability of a positive k was set at a conservative lower level (0.515) derived from simulated populations of constant size (Reich et al. 1999). The interlocus g-test statistics were compared with fifth-percentile values (given by Reich et al. 1999; Table 1) that were appropriate for the number of loci and the sample sizes.
Mitochondrial DNA analysis
Nucleotide diversity (Θπ) and haplotype diversity (h) were estimated using DnaSP version 4.50.3 (Rozas et al. 2003). The Ewens–Watterson homozygosity test for fit of the observed haplotype distribution to the distribution expected under an infinite-alleles model for a population at mutation–drift equilibrium (Ewens 1972; Watterson 1978) was conducted using ARLEQUIN with 10,000 permutations. Observed haplotype numbers (AO), expected haplotype numbers (AE), Tajima's (1989) test of neutrality (DT statistic), and Fu's (1997) test of neutrality (FS statistic) were evaluated under the infinite-sites model with ARLEQUIN, each with 10,000 data permutations to test for significance.
Results
Nuclear Loci
The HO among microsatellite loci ranged between 0.341 and 0.961 (Table 2). No linkage disequilibrium was found for any of the locus pairs. In general, locus–sample combinations conformed to HWE, with only 1 of 72 combinations showing a significant heterozygote deficit. However, global tests of loci over samples revealed heterozygote deficits for Oel42, Pmo367, and Pmo69 (P < 0.001 for all tests). Analysis with Micro-Checker indicated that null alleles were evident in two or three samples for each locus, as well as over all samples combined. Attempts to correct genotype frequencies for null alleles using several correction methods in Micro-Checker or the “including null alleles” (INA) method (Chapuis and Estoup 2007) implemented in FreeNA did not fully resolve these deficits, perhaps due to an inability to discriminate between “true” nulls and PCR amplification failures of one (e.g., through upper allele dropout) or both alleles. Eliminating these three loci resulted in a data set of six loci exhibiting six heterozygote deficits in five of the eight samples, but no significant deficits were found in single-locus tests over individual samples or over all samples for all loci (U-test: P = 0.319).
In exact tests of differentiation between sample pairs, 14 of 28 genic differentiation tests and 11 of 28 genotypic differentiation tests were significant before sequential Bonferroni adjustment, but only three genic and three genotypic differentiation tests remained significant (P < 0.05) after Bonferroni correction for multiple tests (Table 3). All three involved comparison of the 2006 sample from Seguam Pass with other locations. Comparisons between temporal replicate samples from the same locations in Seguam Pass and Stalemate Bank were not significant. Log-likelihood tests indicated heterogeneity over all samples for four loci (Pmo70, Pmo152, Pmo268, and Oel42) and over all loci combined (P = 0.0001).
Despite evidence for patchy variation among locations, the overall level of genetic differentiation between samples was low and nonsignificant. The FST values among sample pairs ranged from 0.0001 to 0.0028, with a mean of 0.0001 for 21 positive estimates. Only three multilocus estimates of FST (ranging from 0.0006 to 0.0021) in 28 pairwise sample comparisons had 95% confidence intervals that did not encompass zero. Similarly, standardized G′ST estimates were quite low, ranging from 0.0000 to 0.0635 and averaging 0.0000 over all loci. The AMOVA results indicated a lack of regional population substructure in a comparison between the Sea of Japan and the Aleutian Islands samples (P = 0.327) and a lack of temporal variability between pooled replicate samples within the latter region (P = 0.672). No significant isolation by distance was detected either across the entire sample range (Figure 2; Spearman's rank correlation: P = 0.824) or among the four Aleutian Islands samples (P = 0.445). Results from STRUCTURE did not infer any geographical subdivision of Atka mackerel across the sampled range. The highest estimated posterior likelihood value for K (loge[−30,932.6]) was obtained for one simulated population, but this value and its variance were not appreciably different from results obtained for specified K-values of 2–6 clusters.
Tests for demographic expansion using microsatellite loci were equivocal. The BOTTLENECK program did not detect departures from neutral allele frequency distributions towards either heterozygosity deficiencies or excesses under a SSM model (P = 0.323) or a TPM model with 5% and 10% multistep mutations (P = 0.156 and 0.119, respectively). In the k- and g-tests, all six microsatellite loci showed negative k-values as expected under a model of recent population growth (P = 0.013). However, the intralocus g-ratio was large (2.142) and nonsignificant.
Mitochondrial DNA
In stark contrast to the large microsatellite and allozyme gene diversities, mtDNA CR sequences showed exceptionally low levels of diversity. Only three haplotypes—including two singletons, each separated by a single substitution from the common haplotype—were detected among 119 individuals from three locations nearly spanning the species' geographic range. These haplotypes yielded very low estimates of Θπ (0.00007) and h (0.033). In addition, mtDNA sequences for a 379-bp fragment of 16S (small subunit) ribosomal RNA (n = 12) and a 799-bp fragment of cytochrome b (n = 10) were monomorphic (data not shown).
Significant departures from neutrality appeared in several tests of mtDNA variability. Over all samples, the AO was 3 and was significantly greater than AE (1.178; P = 0.008). This led to an observed “homozygosity” (0.967) that significantly exceeded expected homozygosity (0.679; P < 0.001). A shift toward low-frequency alleles was also indicated by Tajima's DT (−1.355; P < 0.05) and Fu's FS (−4.294; P < 0.001).
Discussion
Analyses of two molecular marker types showed virtually no genetic structure among populations of Atka mackerel over a substantial part of the species' range. The striking discordance in diversity between nuclear DNA and mtDNA suggests differential genomic responses to historical demographic and evolutionary processes that influence interpretations of contemporary genetic population structure and the phylogeography of this species.
Contemporary Genetic Population Structure
The observed homogeneity for microsatellites among our samples and for allozymes in Lowe et al. (1998) was due to allele frequency homogeneity and not to the analysis of loci with low levels of polymorphism. The distributions of mtDNA haplotypes also did not reveal any geographical structure; however, in this case the power to detect population differences was extremely low because of the very low level of polymorphism. The average H for allozymes (0.137; Lowe et al. 1998) was more than twice that reported by Ward et al. (1994) for 57 marine fish species, and the average H for microsatellites in this study (0.797) was essentially identical to that reported for 12 marine species at 66 loci (H = 0.79; DeWoody and Avise 2000). The low levels of FST and the lack of significant geographic patterns in genetic differentiation of microsatellites (Figure 2) suggested that high levels of gene flow are maintained throughout contiguous portions of the range.
Dispersal ability is inversely related to population genetic differentiation estimated by allozymes and mtDNA in 333 vertebrate species (Bohonak 1999). Adult Atka mackerel do not appear to be highly migratory (McDermott et al. 2005), but the highly protracted spawning–brooding reproductive cycle (up to 6 months) and planktonic larval phase (up to 1 year), combined with vigorous current systems in the Gulf of Alaska and Bering Sea, suggest that long-distance dispersal is likely (Gorbunova 1962; Kendall and Dunn 1985) to impede genetic divergence. The Alaskan populations examined in the present study inhabit two major current systems that are partially isolated by the Alaska Peninsula and Aleutian Islands chain. The inshore Alaska Coastal Current (ACC) flows over a wide continental shelf in the Gulf of Alaska and is influenced by numerous shelf canyons and a complex shoreline. It narrows along the southern edge of the Alaska Peninsula to feed the fast-moving Alaska Stream (Stabeno et al. 2004). A portion of the ACC filters through passes in the Aleutian Islands and drives the eastward-flowing Aleutian North Slope Current and the northwestern Bering Slope Current along the broad, shallow continental shelf in the southeastern Bering Sea (Stabeno et al. 2001). While some of these currents may facilitate larval retention, especially in nearshore areas with complex circulation features (Ladd et al. 2005), larvae entrained in offshore currents are potentially transported along the Aleutian Islands archipelago and into the Bering Sea. Ichthyoplankton surveys have documented Atka mackerel larvae in the Bering Sea (Matarese et al. 2003); these fish were probably spawned in the major nesting areas along the Aleutian Islands, reinforcing the notion that the extent of larval dispersal via fast-moving current systems in both the Gulf of Alaska and the Aleutian Islands archipelago is sufficient to create the broad-scale genetic homogeneity of allozyme and microsatellite loci observed in these regions.
Demography, Natural Selection, and Contrasting Nuclear and Mitochondrial Diversities
A second factor that may have contributed to apparent genetic panmixia in Atka mackerel involves changes in rangewide distributions in response to climate shifts in the North Pacific. The mode and timing of range shifts are relevant for understanding the present-day genetic population structure of many marine species in the North Pacific. In contemporary times (i.e., the period for which oceanographic and biological data exist), this region and adjoining seas have experienced decadal ocean climate shifts, which alter processes influencing primary and secondary productivity (Shuntov et al. 1996; Stabeno et al. 2001) and hence affect abundances of fish at intermediate and high trophic levels (Vasil'kov and Glebova 1984; Francis et al. 1998; Napp et al. 2000; Hunt et al. 2008), such as the walleye pollock Theragra chalcogramma (Stepanenko 1997; Shima et al. 2001). Abundances of several marine fishes have varied by at least an order of magnitude over the past several decades (Spencer and Collie 1997), often in response to climate change. If recent colonization events from genetically homogeneous refuge populations have produced the appearance of broad-scale genetic panmixia, estimates of migration from molecular markers may greatly overestimate the magnitudes of gene flow and the spatial extents of demographically distinctive populations.
More distant historical events associated with ice age climate cycles can also influence the contemporary distributions of genetic population markers. These influences can be inferred by examining the structures of microsatellite trees and mtDNA genealogies for expected shifts in topology as populations mature over long periods (Nei et al. 1975; Maruyama and Fuerst 1984). For example, genetic diversity is lost in small colonizing populations or in populations experiencing a bottleneck. Even in marine species with large population sizes, population declines are expected to lead to the loss of rare alleles despite the fact that H is not diminished (Ryman 1994). If populations decline substantially or if they remain small, H is also expected to decline (Nei et al. 1975). As diploid populations rebound, new mutations accumulate and produce excesses of low-frequency alleles (heterozygosity deficiency) for about 0.4Ne generations (Maruyama and Fuerst 1984). Large populations experience a proportionally greater loss of alleles after a bottleneck (Ryman et al. 1995) and produce a greater excess of low-frequency alleles during recovery (Maruyama and Fuerst 1984). The mutation frequency spectrum asymptotically approaches mutation–drift equilibrium after about 2Ne generations as low-frequency mutations drift to intermediate frequencies.
The evidence for population contraction–expansion as the sole demographic force driving the discordance between nuclear and mitogenomic diversities in this study appears to be equivocal. Intralocus k-tests were significant for each microsatellite locus, inferring that coalescent trees from six independent loci supported a scenario of population expansion for Atka mackerel. However, interlocus g-tests were not significant, perhaps due to mutation rate differences among loci producing different tree topologies, and results from BOTTLENECK did not reveal significant excesses in HE associated with a recent reduction in population size. Our estimates of h (0.033) and Θπ (0.00007) for the mtDNA CR are among the lowest values on record for a marine organism (Orbacz and Gaffney 2000; Garber et al. 2005; Hoelzel et al. 2006; Lu et al. 2006; Carlsson et al. 2007), suggesting the occurrence of a population bottleneck or founder event that greatly eroded mitogenomic diversity but that left nuclear genomic diversities largely unaffected. There is no obvious biological cause (e.g., a highly skewed sex ratio) that would account for this discrepancy; a bottleneck in the mtDNA genealogy would not be affected by the polygynandrous mating system in Atka mackerel, since it is only dependent upon the effective number of females in the population.
An alternative explanation to a bottleneck event for the extremely low levels of mtDNA diversity in Atka mackerel may be selective effects: either a selective sweep (the near fixation of one haplotype due to a fitness advantage) or background selection that keeps slightly deleterious mutations from reaching high frequencies. This latter scenario predicts an increase in diversity with population size (Charlesworth et al. 1995), which was not observed in our data. An analysis of mtDNA sequence variation in 1,638 species representing eight animal taxa (Bazin et al. 2006) supported the hypothesis of selective sweeps resulting in recurrent fixation of haplotypes and loss of variation at linked loci (“genetic draft”; Gillespie 2000). Positive selection can produce similar reductions in diversity and departures from neutrality as a population expansion, and the same star-like haplotype genealogy expected for neutral genes in a growing population can appear as low-frequency haplotypes accumulate around a selected haplotype. While this mode of selection appears to account for patterns of mtDNA variation in some species (e.g., Schizas et al. 2001; Ballard and Whitlock 2004; Ballard and Rand 2005; Zink 2005) and has been inferred for some fishes (Árnason 2004; Grant et al. 2006; Vigliola et al. 2007), the accrual of slightly deleterious mutations in whole-mitogenomic comparisons of several gadid species indicated a nearly neutral pattern of evolution (Marshall et al. 2008).
Demographic expansions and selective sweeps are not mutually exclusive phenomena, and if both processes have occurred, they may reflect population responses to massive ocean climate shifts in the North Pacific during Pleistocene ice ages. Over tens of thousands of years, Northern Hemisphere terrestrial glaciations produced drops in sea level that exposed the shallow continental shelf of the southeastern Bering Sea (Mann and Hamilton 1995; Ager 2003). Populations of Atka mackerel, if present, were probably driven to extinction in the central North Pacific by coastal glaciers or were displaced into southern refugia during periods of cooling (Barrie and Conway 1999; Marko 2004). The most likely postglacial colonization route involved a northeastward expansion of Atka mackerel to the Aleutian Islands archipelago, the contemporary center of its distribution (Logerwell et al. 2005; Lauth et al. 2007ba), from the boreal western Pacific Ocean after glacial retreat in the Holocene. The islands and associated submerged Aleutian platform formed the western tip of the Cordilleran ice sheet during the last glacial maximum (LGM) in the Pleistocene (∼18,000 years before present), although its exact westernmost terminus is unknown (Thorson and Hamilton 1986). The prevailing climate regime during that time rendered the region uninhabitable for Atka mackerel both in temperature and lower trophic level productivity. Present-day temperature in the Aleutian passes ranges from 3.6°C to 8.2°C (Stabeno et al. 2005), and the observed range at inshore spawning sites (3.9–10.5°C; Lauth et al. 2007ba) indicates a preference for warmer waters. Atka mackerel are highly planktivorous, primarily consuming euphausiids and calanoid copepods (Yang 1999), and severe disruption of oceanic circulation patterns at the LGM greatly limited primary productivity in the North Pacific Ocean and Bering Sea (Sancetta et al. 1984).
Paleoclimatic data and limited archaeological evidence suggest that Atka mackerel colonized the Aleutian Islands archipelago between 12,000 and 6,000 years before present. After a thermal maximum at the beginning of the Holocene (approximately 12,000–8,000 years before present), the Aleutian low-pressure cell, which dominates wind patterns in the northeast Pacific, migrated to its present mean position, as did the North Pacific Current, a primary contributor to the Alaska Stream and ACC (Sabin and Pisias 1996). This likely created suitable temperature and productivity regimes for colonization by Atka mackerel before their first documented presence in the Aleutian Islands archipelago. Midden remains from Shemya Island, Alaska, confirm the presence of Atka mackerel at the westernmost end of the Aleutian Islands archipelago at 3,000 years before present (Orchard 1998; S. Crockford, Pacific Identifications, Inc., personal communication). Two middens on Rat Island (western Aleutian Islands) yielded different hexagrammid species compositions for two different age estimates. Identifiable hexagrammid remains in the earlier midden (approximately 1,240 years before present) were dominated (88%) by the rock greenling Hexagrammos lagocephalus and kelp greenling H. decagrammus, whereas hexagrammid remains from the more recent site (about 720 years before present) consisted mostly (54%) of Atka mackerel (M. Partlow, Central Washington University, Ellensburg, personal communication). This observation is consistent with a putative eastward colonization path through the Aleutian Islands from a refuge near the species' origin. The genus Pleurogrammus is postulated to have arisen in the western boreal Pacific (Shinohara 1994), eventually giving rise to the Atka mackerel and arabesque greenling P. azonus, which are currently sympatric in the southern part of the Sea of Okhotsk, the Kuril Islands, and the Pacific coast of Hokkaido, Japan. While a comprehensive view of Atka mackerel colonization patterns across the Aleutian Islands is lacking, the presence of this species in the Aleut paleodiet implies that it was locally abundant dating back at least 3,000 years.
Results from this study give some insight into the potential mode of postglacial colonization of the central North Pacific. Previously described patterns of diversity resulting from “pioneer” or “phalanx” colonization (Hewitt 1996) are not apparent in Atka mackerel. Latitudinal gradients in genetic diversity resulting from pioneer colonizations by a few individuals, as seen in high-latitude lake fishes (Bernatchez and Wilson 1998), were not observed. The high levels of polymorphism observed in microsatellite loci are consistent with a large-scale phalanx colonization scenario, but the reduced levels of mtDNA diversity reject this model unless selection is invoked. Conformation to a “center–periphery” model in which marginal populations are expected to be less genetically diverse is not evident our data, although additional samples collected nearer to the current range periphery in the northwestern and northeastern Pacific Ocean would be necessary to more rigorously test this hypothesis. The large spatial extent of low mtDNA diversity observed in Atka mackerel could potentially have resulted from population bottlenecks within refugia, colonizations by only a few individuals, a virtually rangewide selective sweep on mtDNA, or some combination of these processes. Postglacial colonization of the central North Pacific is likely to have occurred with Holocene warming of the region, and newly founded populations have subsequently grown to large contemporary abundances in the Aleutian Islands—a conclusion supported by mtDNA haplotype departures from neutrality as measured by Tajima's DT and Fu's FS.
This paradox of contrasting nuclear and mtDNA diversities supports the conclusions from theoretical (Ryman et al. 1995) and empirical (Smith et al. 1991; Pichler and Baker 2000; Hauser et al. 2002) studies showing that apparently large populations of marine species are not immune to the loss of genetic diversity, although the processes responsible for these losses may be difficult to decipher. A scenario in which Ne falls below thresholds that purge nearly all variation in mtDNA across most of the species' geographic range but does not fall low enough or for a sufficient duration to affect diversities of nuclear markers seems unlikely without invoking natural selection. If both neutral and selective forces have been important in producing discordant patterns of gene diversity, they are both probably tied to historical ocean climate changes across the region.
Acknowledgments
We thank G. Jensen, B. Flerx, N. Roberson, and J. Orr for sample collection. This research is Contribution EcoFOCI-R745-RAAO to National Oceanic and Atmospheric Administration (NOAA) Fisheries Oceanography Coordinated Investigations. Susanne McDermott and Kim Rand provided helpful comments on earlier drafts of this manuscript. Terry Josey kindly provided the drawing of Atka mackerel displayed in Figure 1. This study was funded in part by the Resource Assessment and Conservation Engineering Division and the Resource Ecology and Fisheries Management Division of the NOAA Alaska Fisheries Science Center. The views expressed herein are those of the author(s) and do not necessarily reflect the views of NOAA or any of its subagencies. Use of trade names does not imply endorsement by the U.S. Government.
References
Figure 2.
Linearized genetic differentiation index FST (calculated as FST/[1 − FST]) versus the natural logarithm of geographic distance between Atka mackerel samples collected in the Sea of Japan and in Alaska (see Table 1 for description of sample locations).
Table 1.
Sample location, collection year, abbreviation, geographic coordinates, and sample sizes for microsatellites (msat) and mitochondrial DNA (mtDNA) of Atka mackerel.
Table 2.
Sample numbers (n), number of alleles (A), allelic richness (AR), expected heterozygosity (HE), observed heterozygosity (HO), and inbreeding coefficient (FIS) values of microsatellite loci in Atka mackerel for each sample and over all samples. For FIS values, significance before (italic) and after (bold) correction for multiple tests is indicated. Sample abbreviations are defined in Table 1.
Table 3.
Probability (P) values for exact tests of allelic variation (below diagonal) and genotypic variation over all microsatellite loci in Atka mackerel. Significant comparisons before (italics and bold) and after (bold only) sequential Bonferroni correction for multiple tests are indicated (initial α = 0.0018). Sample abbreviations are defined in Table 1.
