Mitochondrial Cytochrome b DNA Variation in the High-Fecundity Atlantic Cod: Trans-Atlantic Clines and Shallow Gene Genealogy

Corresponding author: Einar Árnason, Sturlugata 7, 101 Reykjavik, Iceland., einararn{at}hi.is (E-mail)

Communicating editor: M. AGUADÉ

	ABSTRACT

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

An analysis of sequence variation of 250 bp of the mitochondrial cytochrome b gene of 1278 Atlantic cod Gadus morhua ranging from Newfoundland to the Baltic shows four high-frequency (>8%) haplotypes and a number of rare and singleton haplotypes. Variation is primarily synonymous mutations. Natural selection acting directly on these variants is either absent or very weak. Common haplotypes show regular trans-Atlantic clines in frequencies and each of them reaches its highest frequency in a particular country. A shallow multifurcating constellation gene genealogy implies young age and recent turnover of polymorphism. Haplotypes characterizing populations at opposite ends of the geographic distribution in Newfoundland and the Baltic are mutationally closest together. The haplotypes are young and have risen rapidly in frequency. Observed differentiation among countries is due primarily to clinal variation. Hypotheses of historical isolation and polymorphisms balanced by local selection and gene flow are unlikely. Instead the results are explained by demic selection of mitochondria carried by highly fit females winning reproductive sweepstakes. By inference the Atlantic cod, a very high-fecundity vertebrate, is characterized by a high variance of offspring number and strong natural selection that leads to very low effective to actual population sizes.

SHALLOW intraspecific gene genealogies characterize many marine fishes although sister taxa often show considerable genetic divergence (GRANT and BOWEN 1998 ; GRAVES 1998 ). This indicates a recent age of extant populations and results from relatively small effective population sizes (N_e), either historical or contemporary, or both. A debate is ongoing about whether historical or contemporary forces are mainly responsible. A historical hypothesis predicts temporally stable spatial divergence among populations whereas an alternative contemporary hypothesis emphasizes temporal aspects of variation. Demographic population fluctuations reduce effective relative to actual population sizes (VUCETICH et al. 1997 ), as do skewed sex ratios. Another important factor reducing effective population sizes is a greater than Poisson variance in offspring number.

Many or most marine organisms differ from the usual model organisms by high fecundities and high mortalities of early life stages or type III survivorship. High fecundities of males may also imply sperm competition. With broadcast spawning and external fertilization, dilution may under some conditions lead to sperm limitation (LEVITAN and PETERSEN 1995 ; YUND 2000 ) and hence competition among eggs for sperm. From these facts results the potential for a high variance in the number of surviving offspring from the parents of the current generation (HEDGECOCK 1994 ) or a potential for a large variance in individual fitness. The variance in net fitness of individuals can be looked at as the variance in the realized number of offspring an individual contributes to the next generation. Fitness of an individual arises from a unique interaction of this individual's genes and its environment during its life from birth to death. However, the genotype of each individual of an outcrossing sexual species is unique and, therefore, the total genotype of an individual cannot be a target for natural selection because individuals do not have heritability (ARNASON and BARKER 2000 ). Instead, specific traits or genes contributing to these traits are the targets of selection such as traits for selection of mates and right spawning conditions. However, other genes in the genome will experience the effects of the fitness differences at target loci as a higher rate of random genetic drift because selection will generate a higher than Poisson variance in offspring numbers.

Even if a locus is not linked to another gene that is a target of selection and the variation remaining at this locus is neutral to selection the locus in question will nevertheless experience the effects of variance in net fitness as a decrease in effective population size and higher rate of random genetic drift. The variance effective number of a bisexual population as a function of the actual number of breeding individuals, N_a, and the mean, , and variance, ²_k, in number of contributed offspring is N_e = (4N_a – 4)/( + ²_k) (WRIGHT 1969 ). Under long-term stability of population size = 2 and with Poisson distribution of offspring numbers the variance equals the mean and, therefore, the ratio of effective numbers to actual numbers of breeders would be 1, N_e/N_a 1. At another extreme with a winner-take-all sweepstakes reproduction (WILLIAMS 1975 ; HEDGECOCK 1994 ) if all surviving offspring come from a single individual the mean number of offspring stays the same but the variance becomes ²_k = N_a and therefore N_e 4 and the N_e/N_a ratio becomes a tiny fraction (FRANKHAM 1995 ; TURNER et al. 2002 ). The variance of offspring numbers among high-fecundity organisms in the real world can potentially lie somewhere on the scale between 2 and N_a.

Instead of looking forward in time the coalescent (e.g., NORDBORG 2001 ) takes the present, the data, and traces the genealogy to coalescence. A simple coalescent models a haploid or clonal organisms of N individuals with a binomial distribution having mean and variance of one, a bifurcating genealogy and coalescent time t = N/², where t = 1 refers to N generations ago. Applying this concept to a winner-take-all sweepstakes reproduction ² = N with a multifurcating genealogy that coalesces immediately in the previous generation and t = 1 refers to one generation ago. Real populations lie somewhere on that continuum. The question of where on that continuum a species lies has implications for the nature of selection, adaptation, and population structure in that species. The contrast drawn here between an organism with Poisson variance in offspring number vs. a winner-take-all sweepstakes makes different predictions about the shape of the genealogy. From the shape of the gene genealogy one may therefore draw inverse inference about variance in fitness although a formal direct estimation of the variance is not possible (NORDBORG 2001 ). The coalescent assumes bifurcating genealogies and breaks down under multifurcating genealogies but coalescences for time-inhomogeneous environments are beginning to be studied (MOHLE 2002 ). High-fecundity marine organisms with fluctuating population sizes would fit in with that. They may have effective population sizes orders of magnitude lower than actual population sizes, most likely due to high variance in offspring number (HEDGECOCK et al. 1992 ; TURNER et al. 2002 ).

High-fecundity marine organisms typically reproduce under spatially and temporally varying oceanographic conditions that affect their sexual maturation, choice of mate, external spawning and fertilization success, survival of gametes, zygotes, and larvae, and recruitment (HEDGECOCK 1994 ). Such fitness differences can result in intense fecundity as well as viability selection (HARTL and CLARK 1989 ). Most of these organisms have planktonic larvae, which is conducive to gene exchange over widely separated areas. Some are highly mobile animals able to travel as individuals of reproductive age over large distances, thus adding to the potential for short- and long-distance gene flow over large areas. An essential feature of such dispersal is that "recruits to a local population come from somewhere else" (JOHNSON and BLACK 1984 , p. 1371). Adaptations to local conditions, therefore, are not expected to accumulate over time and changes in genetic composition are due to short-term or even single-generation effects of selection and recruitment (WILLIAMS 1975 ; JOHNSON and BLACK 1984 ). There can nevertheless be genetic patchiness, or structure, in time and space but if reproduction is characterized by sweepstakes that by chance match oceanographic conditions favorable for individual fitness (HEDGECOCK 1994 ) any such structure will be ephemeral. However, there can be some memory in the system that sets a timescale, for example, large distances between local spawning sites or patchy oceanographic conditions suitable for reproduction.

The Atlantic cod, Gadus morhua, is a high-fecundity marine organism with highly mobile adults, external fertilization, and planktonic larvae. The opportunity thus exists in cod for considerable gene flow and also for a high variance in offspring number or a high variance in net fitness. A question arises of how genetic variation is molded in high-fecundity species (TURNER et al. 2002 ) such as cod, at loci that are direct targets of selection as well as at loci influenced by selection at other genes through either linkage or reduction in effective population sizes relative to actual population sizes.

Studies of nuclear DNA variation in Atlantic cod, both microsatellites (BENTZEN et al. 1996 ; RUZZANTE et al. 1996A , RUZZANTE et al. 1996B , RUZZANTE et al. 1997 , RUZZANTE et al. 1998 ; RUZZANTE 1998 ) and restriction fragment length polymorphisms (RFLPs; POGSON et al. 1995 , POGSON et al. 2001 ) and protein-coding genes (POGSON 2001 ), have been interpreted to show isolation by distance on large scales as well as microspatial population differentiation and temporal stability on both annual (JONSDOTTIR et al. 2001 ; BEACHAM et al. 2002 ;) and decadal (RUZZANTE et al. 2001 ) scales. Mitochondrial DNA (mtDNA) studies have not found evidence for significant population structure within cod stocks of single countries (CARR and MARSHALL 1991A ; PEPIN and CARR 1993 ; CARR et al. 1995 ; ARNASON and PALSSON 1996 ; ARNASON et al. 1998 , ARNASON et al. 2000 ; SIGURGISLASON and ARNASON 2003 ), yet they have found a small but significant temporal effect in some instances (ARNASON et al. 1998 , ARNASON et al. 2000 ). This study is an analysis of the entire data set of mtDNA variation of Atlantic cod from the North Atlantic. It is expected that this large study can reveal new patterns in the data and understanding of the forces responsible for reducing effective relative to actual population numbers. Using nested clade analysis it attempts to make inference about historical and contemporary forces (TEMPLETON et al. 1995 ; TEMPLETON 1998 ).

	MATERIALS AND METHODS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Previous articles have described mtDNA variation of Atlantic cod stocks from single countries. Here an extensive analysis is made of the entire mtDNA cytochrome b data set on Atlantic cod.

The data:
The data are a sequence variation of a 250-bp cyt b fragment corresponding to base pairs 14459–14708 included in the JOHANSEN and BAKKE 1996 complete mtDNA sequence. The samples come from various localities throughout the Atlantic, ranging from Newfoundland (CARR and MARSHALL 1991A , CARR and MARSHALL 1991B ; PEPIN and CARR 1993 ; CARR et al. 1995 ), Greenland and Iceland (ARNASON et al. 2000 ), the Faroe Islands (SIGURGISLASON and ARNASON 2003 ), and Norway (ARNASON and PALSSON 1996 ) to the Baltic and White Seas (ARNASON et al. 1998 ).

Molecular analysis:
The methods of molecular analysis followed those already described (ARNASON and PALSSON 1996 ; ARNASON et al. 1998 , ARNASON et al. 2000 ; SIGURGISLASON and ARNASON 2003 ). Briefly, DNA was isolated using methods ranging from genomic DNA isolation to simple proteinase K digestion of a small amount of tissue. Subsequently DNA was amplified by PCR, sequencing reactions were done either manually or by thermal cycling, and electrophoresis also was done either manually or by automatic sequencers as already described. The 250-bp fragment was common to all the studies.

Genealogy:
To facilitate the evaluation of plausible gene genealogies, parsimony trees were constructed using PHYLIP (FELSENSTEIN 2002 ), PAUP (SWOFFORD 2002 ), and TCS (CLEMENT et al. 2000 ), which implements statistical parsimony. To further evaluate which one of the alternative genealogies was most plausible, the frequency of haplotypes was used as an additional criterion. Thus rare and singleton haplotypes were derived from more common haplotypes in accordance with predictions from coalescence theory (POSADA and CRANDALL 2001 ). Using this criterion also in most cases happened to coincide with geographical distribution so that rare haplotypes in conflict were derived from a haplotype from the same population, a criterion in accordance with molecular variance parsimony (EXCOFFIER and SMOUSE 1994 ).

Analysis of molecular variance and population structure:
Population genetic structure was inferred by analysis of molecular variance (AMOVA) using the Arlequin package (SCHNEIDER et al. 2000 ). Variation was partitioned using TAMURA and NEI's (1993) genetic distance between haplotypes and also by considering haplotypes as names, discounting only their genetic differences. Neutrality was tested and various neutral theory statistics of diversity and effective population numbers scaled by mutation rate, or , were estimated with Arlequin. Population differentiation as pairwise F_ST and gene flow as effective number of migrants per subpopulation per generation (SLATKIN 1995 ) were also estimated.

Nested clade analysis:
The GeoDis package (POSADA et al. 2000 ) was used for nested clade analysis (NCA; TEMPLETON et al. 1995 ; TEMPLETON 1998 ) of geographic association of haplotypes and for a geographic distance analysis by permutation. NCA applies the genealogical information for nested statistical tests, identifying whether patterns are due to contemporary forces of isolation by distance or historical events. The method estimates clade dispersion, D_c, which measures the geographical range or dispersion of a particular clade and the nested clade distance, D_n, which measures geographical distribution or displacement of a particular clade relative to other clades nested with it at a particular level of the hierarchy. The method defines interior (I) and tip (T) clades of a genealogy and estimates distances for the I – T difference. Inference on historical or contemporary forces can be drawn with a standardized key (TEMPLETON et al. 1995 ) on the basis of whether distances are significantly small or large.

Historical demography and ratio of effective to actual population numbers:
Scaled effective population size was estimated with the program fluctuate (KUHNER et al. 1995 , KUHNER et al. 1998 ) under a model of no growth (_t = ₀) and under a model of exponential growth (_t = ₀e^gt). Here = 2N_eµ and g is in units of µ^–1. Typical runs involved 10–20 short chains and 2 long chains of at least 6000 and 12,000 steps, respectively, with a 30-step sampling increment. Statistical contour descriptions of maximum log-likelihood obtained under exponential growth were made using polynomial trend surface regression (VENABLES and RIPLEY 2002 , p. 420). Further exploration of parameter space was done, including restarts at or near the local peaks obtained and using the same number of steps and chains. Some restarts were done with 2–4 short chains and 2–4 long chains of 18,000 steps. Historical demography was further explored with a generalized skyline plot, which is a graphical nonparametric model-free estimate of through time embodied in a genealogy (STRIMMER and PYBUS 2001 ). The method assumes a clocklike behavior and a UPGMA tree was made using PHYLIP and used for skyline estimation performed under R with the APE package (http://cran.r-project.org). The UPGMA tree was also used as an initial tree for estimation with fluctuate. The ratio of effective to actual population numbers due to fluctuating population sizes was evaluated using a method of VUCETICH et al. 1997 with fisheries data of Icelandic cod between 1955 and 2002 (ANONYMOUS 2003 ).

To estimate effective size N_e from , mutation or substitution rates are required. A direct estimate of cod mtDNA mutation rate is, with 95% confidence, < 1.5 x 10^–5 (SIGURGISLASON 2003 ; H. SIGURGÍSLASON and E. ÁRNASON, personal communication). A calibration of gadid mtDNA molecular clock using the trans-Arctic interchange (VERMEIJ 1991 ) and molecular data (CARR et al. 1999 ) yields a synonymous substitution rate of = 3.86 x 10^–8/site/year (E. ÁRNASON, unpublished data), which, with a generation time of 4.8 years (E. ÁRNASON, unpublished data), gives a substitution rate of = 1.85 x 10^–7/site/generation. Some results and analysis are made available as supplemental data on the Genetics website at http://www.genetics.org/supplemental/.

	RESULTS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Molecular and genealogical relationships:
A total of 39 sites were found segregating for variation among the 1278 individuals defining 59 observed haplotypes (Fig 1). All but seven mutations were synonymous. Nonsynonymous mutations were observed at site 467, defining haplotypes U and UI; sites 468, 522, and 542, defining haplotypes ZI, F, and L, respectively; and site 695, defining haplotypes PI, MC, and M. To account for these amino acid variants the same site must have been hit more than once. Most amino acid replacements were functionally conservative Isoleucine for Valine or Valine for Isoleucine replacements found as singletons in the sample. Atlantic cod differed from the related Greenland cod, G. ogac, by 15 synonymous substitutions (CARR et al. 1999 ). The apparent excess of nonsynonymous intraspecific mutations to interspecific substitutions (MCDONALD and KREITMAN 1991 ), however, was not significant (Fisher's exact test, P = 0.33; similarly counting sites instead of mutations, 5/34 vs. 0/15, P = 0.31).

View larger version (31K): In this window In a new window Download PPT slide   Figure 1. Segregating sites of 59 haplotypes of a 250-bp fragment of cytochrome b gene found among 1278 cod sampled combined from various localities in the North Atlantic. The site number is from the complete cod mtDNA sequence (JOHANSEN and BAKKE 1996 ) less 14,000. N is the frequency of a haplotype in this combined sample. A dot is identity of sites with the A haplotype reference sequence. A star indicates an amino acid variant.

High-frequency haplotypes (Fig 2) were genetically close to each other with E, G, and C one mutation away from the most frequent haplotype A, and D, derived from C, two steps away. They were all internal to the tree and a number of rare or singleton haplotypes were derived from each high-frequency type, or, alternatively, a number of rare types simultaneously coalesced into a common type. Frequency of common haplotypes and the number of rare haplotypes derived from each were significantly correlated (Spearman's rank correlation, P = 0.03). The high-frequency haplotypes thus defined clades as being "stars" in a "constellation" phylogeny. Four of the less common haplotypes (NI, H, DI, and B) also spawned off rare haplotypes. Only a single internal node in the network was not represented by an individual in the sample.

View larger version (29K): In this window In a new window Download PPT slide   Figure 2. Phylogenetic maximum-parsimony network of 59 haplotypes of a 250-bp fragment of cytochrome b gene found among 1278 cod sampled from various localities in the North Atlantic. Unit length of line is change of 1 bp.

A large number of equally parsimonious 60-step trees were found and considerable homoplasy was apparent (Fig 1). Most conflicts could be resolved by parsimony. The remaining conflicts were resolved by a frequency criterion. Haplotypes LI and LJ provided an example. Both shared a purine transition at site 487 but LI also shared a transition at site 631 with haplotype NI. Therefore LJ could have been derived from LI with a back transition at site 631 or it could have been derived from LJ with an independent gain of the transition at site 631. However, haplotypes G and NI are more frequent than both LI and LJ. Therefore, the genealogy was resolved as shown in Fig 2 by deriving LI from the more common NI and not from the rare haplotype LJ. The frequency criterion also coincided with a criterion of geography because the rare or singleton haplotypes that showed conflicts in most cases were not found in the same population. Applying these additional criteria fully resolved conflicts (Fig 2). The haplotypes were grouped into five clades (Fig 2) that form a basis for nested clade analysis.

The transition/transversion ratio was 29:1 (and see Table S1 at http://www.genetics.org/supplemental/). The two types of relative purine transition rates were similar as were the pyrimidine transition rates. However, the overall relative frequency of purine transition mutation was about twice as high as the frequency of pyrimidine transition (Table S1). In contrast, the interspecific purine:pyrimidine transition ratio for cytochrome b between Atlantic cod and Greenland cod was 3:12 and for Atlantic cod and Walleye pollock, Theragra chalcogramma, 2:13 (CARR et al. 1999 ). In both instances tests of the intra- vs. interspecific ratios were significant (Fisher's exact test, P = 0.02 and P = 0.004, respectively).

The 82 codons had five guanine synonymous sites. Ten mutations were observed at four of the guanine sites, or 2.5 mutations per segregating synonymous guanine site, considerably higher than the overall hit rate of 1.54. There were also significant differences in the synonymous purine:pyrimidine intraspecific variation of Atlantic cod in comparison with interspecific differences in Greenland cod and Walleye pollock (Fisher's exact test, P = 0.04 and P = 0.016, respectively).

There was an apparent heterogeneity of nucleotide mutation frequency among sites. On the basis of the genealogy as drawn (Fig 2), the observed number of hits per site (based on all sites) had a mean and variance of 0.24 and 0.43, respectively. The observed numbers did not fit a Poisson distribution (G = 39.3; P < 10^–6; d.f. = 3) but gave a satisfactory fit to a negative binomial distribution (G = 3.1; P = 0.21; d.f. = 2; Table S2 at http://www.genetics.org/supplemental/). On the basis of this and because of a difference in the purine and the pyrimidine mutation frequencies, the model of nucleotide mutation discussed by TAMURA and NEI 1993 was considered appropriate for AMOVA with a gamma distribution parameter = 0.30 (obtained by method of moments; RICE 1995 , p. 278).

Trans-Atlantic clines in haplotype frequencies:
The high-frequency haplotypes were found in all countries except in the small sample from the White Sea. The C haplotype similarly was found in all countries, but NI, which also exceeded 1% frequency, was found in all countries except Newfoundland. There were highly regular trans-Atlantic clines in the frequencies of the high-frequency haplotypes or stars of the constellation (Fig 3). Geographically intermediate localities as a rule had intermediate haplotype frequencies. Haplotype A, overall the most common, reached its highest frequency in Newfoundland, where it was the dominant haplotype, and decreased in a regular manner toward the east. Haplotype E, overall the second most common, reached its highest frequency in the Baltic. It was also the most frequent haplotype in the Baltic, although it was not a dominant haplotype in the manner that haplotype A was in Newfoundland. Haplotype E decreased in frequency going westward away from the Baltic. In a similar fashion haplotype G reached its highest frequency in Norway and haplotypes D and C reached theirs in Iceland and all decreased in frequency in a regular fashion going away from these localities. The clines of the different high-frequency haplotypes were statistically independent. For example, considering the four high-frequency haplotypes and binning the rest over six localities, there were 20 d.f. and thus there was ample room for independent variation. The regularity of clines thus was not due to statistical confounding. The clines for the major clades (Fig 2) were very similar to clines of frequencies of high-frequency haplotypes (Fig 3).

View larger version (16K): In this window In a new window Download PPT slide   Figure 3. Clines in frequency of high-frequency haplotypes on an ordinal scale from Newfoundland in the west to the Baltic in the east. mtDNA cytochrome b sequence variation among 1278 Atlantic cod is shown.

Tests of neutrality:
With the exception of Newfoundland, both haplotype and nucleotide diversities were high in all countries and in the overall sample (Table S3 at http://www.genetics.org/supplemental/). On that account countries except Newfoundland fall into a category of high and (GRANT and BOWEN 1998 ). The reason, however, is not deep divergence but relatively high and even frequencies of many haplotypes. The high and dominant frequency of the A haplotype in Newfoundland was mainly responsible for the low diversities there because frequencies of haplotypes influence these statistics. The differences among countries were smaller in the estimates of based on number of haplotypes and number of segregating sites, which are governed by sample size and not by frequencies of haplotypes. TAJIMA's (1989) D was highly significant and negative for Newfoundland and nonsignificant or at borderline for other countries. Iceland had a negative D at a borderline of significance most likely due to the large sample size in Iceland. A large sample picks up rare haplotypes, which increases _S more than , resulting in a negative D. This effect of sample size of picking up rare haplotypes also was seen in the _k values. FU's (1997) F_S statistic, which is sensitive to the presence of rare mutations relative to expectation based on nucleotide diversity, was highly significant for Newfoundland, Iceland, Greenland, and the Faroe Islands. However, Norway and the Baltic Sea, which were farthest away from Newfoundland, which had a significant D, were not significant for F_S.

AMOVA and partitioning of variation:
Areas within countries explained none of the variation by AMOVA (Table 1, bottom). The resulting correlation of random haplotypes within areas of countries relative to that of random haplotypes drawn from the country to which those areas belonged (_SC) had a negative sign and was interpreted as nil. The correlation of random haplotypes within areas of countries relative to random haplotypes drawn from the whole (_ST) and the correlation of random haplotypes within countries relative to random haplotypes from the whole (_CT) were very similar and both were positive and significant. Since the effect of areas within countries was very small it was ignored and the AMOVA was done on differences among countries only (Table 1, top and middle). Using TAMURA and NEI's (1993) genetic distances among haplotypes, differences among countries explained 8% of the variation. However, discounting genetic differences among haplotypes and treating them as names only, countries explained almost twice as much or 14% of the variation. The among-countries variance was almost the same and thus was not affected by discounting genetic differences. However, the within-countries component of the variance was reduced one-half by discounting the genetic differences among haplotypes, resulting in a higher percentage explained. Thus the differentiation among countries was due primarily to the frequency differences among haplotypes, particularly the clinal variation.

Pairwise population differences and gene flow:
Newfoundland and the Baltic were mainly responsible for differences among countries in terms of estimated F_ST and gene flow (Table 2). Both countries showed significant differences from all the others. However, regular clines also were evident in these statistics, probably due to the haplotype-frequency clines. There was an order of magnitude step going from Newfoundland and the Baltic to the four central countries, Greenland, Iceland, the Faroe Islands, and Norway. The estimated gene flow was greater than or equal to one, which is a cutoff sufficient to prevent differentiation by random drift alone under the island model. Although Greenland and Norway showed weak differences (P = 0.044 ± 0.003) there was considerable gene flow among the central countries. Pairwise F_ST had negative signs between Iceland and the Faroe Islands and the Faroe Islands and Norway, which were interpreted as nil, resulting in infinite N_em.

Genealogy, geography, and nested clade analysis:
Applying a definition of clades (Fig 2), the haplotype distributions were formally tested, using a nested clade contingency analysis (Table S4 at http://www.genetics.org/supplemental/). The rare tip clades (Fig 2) showed either no geographic or no genetic variation and thus were not tested (TEMPLETON 1998 ). Of the single-step clades only D showed significant geographic association under contingency tests. Combining C and D in a two-step clade, 2-DC, however, showed no significance. The two-step clade 2-AEG and the total cladogram showed highly significant geographic association.

A nested clade geographic distance analysis (Table 3) was made to further explore and determine if the patterns resulted from contemporary or historical forces. Haplotypes or zero-step clades were tested within one-step clades and only the high-frequency haplotypes and haplotypes showing significant distances are presented in the table. The high-frequency haplotypes were interior nodes within each one-step clade. A caveat is introduced that numerous tests were performed on haplotypes within one-step clades and only a few showed significance, mostly at a level of 0.05 > P > 0.01. After a Bonferroni adjustment (RICE 1995 ) only two remained significant at the zero-step level. Four rare haplotypes within the A clade, N, AJ, RI, and X, were found in a single country and showed significantly small D_c distances. They were too frequent within their respective countries and would be expected to be more dispersed under the null hypothesis. Haplotype MI, found in one individual from Greenland, three from Iceland, and three from the Faroe Islands, similarly was restricted in distribution relative to its frequency as shown by a significantly small D_n distance. Haplotype B, however, found in Newfoundland and Norway, showed a significantly large D_n distance. The dispersion distance D_c was significantly large for the interior vs. tip comparison. With the caveat above, the standardized inference was that haplotypes within the A clade showed both restricted gene flow and long-distance dispersal. Within the E clade the interior E haplotype was overdispersed (large D_c) and also overdisplaced relative to other haplotypes within the clade (large D_n). The inference was restricted gene flow and isolation by distance. The null hypothesis was not rejected within the G clade. Within the D clade haplotype J had a significantly large displacement (D_n), which, coupled with a significantly short I –T distance, led to an inference of range expansion and long-distance colonization. Similar range expansion was inferred within the C clade.

At the next level within the two-step clades the null hypothesis was not rejected for clades D and C. Within the AEG clade, however, clade A showed a significantly large distance displacement (D_n) from the center of the AEG clade. The E and G clades, however, both showed significantly small dispersion (D_c) and the G clade also showed a significantly small displacement relative to the two other clades within AEG. With the A clade forming an interior node of AEG, the I – T distances were both significantly large. Similarly within the total cladogram both distances were significantly large for AEG but significantly small for the DC clade and for the I – T comparison. The inference was restricted gene flow and isolation by distance for both AEG and the total cladogram.

Effective population size:
An estimate of = 0.0046 had the highest log-likelihood under the model of no growth (Fig 4A). The highest log-likelihood minus 1/2²_1,0.05 is an 95% confidence limit (dotted line). The log likelihoods for the lowest and highest that almost reached or exceeded this limit were obtained for of 0.0032 and 0.0092, respectively, thus defining a range of values with some support. The generalized skyline plot nonparametric estimate of through time (STRIMMER and PYBUS 2001 ) gave similar results to the no-growth model (Fig 4B). The skyline showed a modest growth in the first half of the shallow genealogy followed by fluctuations around some median value and an increase in the most recent time resulting, however, in a limited population size ( = 0.036). Skyline plots of representative best trees obtained with fluctuate, however, showed the opposite effect of a sharp decline in the present (data not shown). The skyline median = 0.0044 was in good accordance with scaled population size estimates based on the mean number of pairwise differences, = 0.0066, but was less than _k and _S, which are sensitive to an excess of rare alleles (Table S3). The estimates of and growth rate g under the model of exponential growth can be grouped into three clusters (Fig 4C). First, with the lowest log-likelihood <2, a cluster of results with high population size and high positive growth is shown (median = 0.16 and = 3066). Second, a cluster with intermediate log-likelihoods, with a population size similar to the no-growth model and slightly positive and slightly negative growth, is shown (median = 0.0033 and = – 367). The log-likelihoods for this cluster were similar but slightly lower than the log-likelihoods under the no-growth model. Finally, and with the highest log-likelihoods between 8 and 45, a cluster of very small and negative growth is shown (median = 0.0006 and = – 5771). A further exploration of this region of parameter space (Fig 4D) showed that the first solution (high population size and positive growth) was unstable. With small perturbations from peaks obtained in first round of estimation the fluctuate estimator yielded results in either the second or the third cluster. In turn, the second cluster could also, albeit with greater difficulty, be perturbed to yield still smaller population sizes and negative growth. Thus the exponential model provided a better fit than the no-growth model only for negative growth.

View larger version (33K): In this window In a new window Download PPT slide   Figure 4. Historical demography of effective population size and growth rate based on mtDNA cytochrome b sequence variation among 1278 Atlantic cod. (A) Log-likelihood on obtained in various runs of a no-growth model (t = 0) with fluctuate. Dotted horizontal line is drawn 1/221,0.05 below the maximum log-likelihood obtained. Dashed vertical line is drawn at median. (B) Generalized skyline plot of log through time from present to past obtained with a nonparametric estimate. (C) Observed parameter combinations of growth rate and obtained under various runs of an exponential growth model (t = 0egt) with fluctuate. Contours are a statistical trend surface description of log-likelihoods of peaks. (D) Representative initial (open circles) and final (solid circles) parameter combination estimates ( and ) obtained under an exploration of the log-likelihood landscape of the exponential growth model in C.

Under mutation-drift steady state the estimates are effective population numbers, N_e, scaled by mutation rate per site per generation. Using the above-mentioned substitution rate of 1.85 x 10^–7 per site per generation, the median N_e would be 432,432, 8919, and 1622 for the three clusters of peaks in parameter space, respectively. The first one was unstable. Similarly they would be 12,432 and 11,892 under the no-growth and skyline estimations, respectively. Using a direct estimate of mutation rate the effective sizes would be smaller still.

The ratio of effective to actual population numbers based on ecological data on population fluctuations (VUCETICH et al. 1997 ) in the Icelandic stock (ANONYMOUS 2003 ) was estimated as N_e/N_a = 0.82. The largest yearly catch of Atlantic cod was 2.7 x 10⁶ tons in 1956 (JONSSON 1992 ). From this the actual population size is not <10⁹ and probably one or two orders of magnitude larger. Thus, on the basis of the above genetic estimates the N_e/N_a ratio is very small (<10^–5–10^–6).

	DISCUSSION

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Neutrality and selection:
The variation observed was mostly synonymous mutations. From the nature of the genetic code it follows that about one-fourth of single-base mutations are to synonymous codons (KING and JUKES 1969 ). Thus the variation observed in cod is likely to be what is left over after purifying selection has taken its toll.

Commonly, more amino acid variants are segregating within species than are found in interspecific comparisons for various genes and organisms, indicating that some weak purifying selection remains (BALLARD and KREITMAN 1994 ; Rand 2001 ). The presence of rare amino acid variants in Atlantic cod relative to their absence in interspecific comparisons, even if not statistically significant, may also be indicative of weak selection. The higher purine-to-pyrimidine transition mutation ratio within cod is due primarily to a high g a bias (and see SIGURGISLASON and ARNASON 2003 ). A possible source of the mutational bias is deamination of cytosine on the opposite strand, perhaps caused by oxidizing radicals from respiration, and subsequent mutation during DNA replication. The interspecific comparison, however, indicates that there is weak purifying selection also at synonymous sites. The TAJIMA 1989 and FU 1997 test results, which were significant in some instances, are compatible with this. These tests are sensitive to excess of rare alleles.

Gene genealogy:
The gene genealogy is shallow as observed in many marine fish (GRANT and BOWEN 1998 ), with high diversity resulting not from deep divergence but from high and relatively even frequencies of many haplotypes, which are internal nodes of the genealogy. The high-frequency types are stars in a constellation genealogy. Almost all internal nodes exist as individuals in the sample. The external nodes are rare, many found only as singletons. This contrasts with mtDNA control region variation in humans, which has several deep branches with many and even most haplotypes found only as external nodes and many or even most internal nodes not found as individuals in the samples (see, for example, minimum spanning networks in BANDELT et al. 1995 ; WATSON et al. 1997 ; MACAULAY et al. 1999 ). The implication is that the cod mtDNA genealogy is young compared to the human and yet it has a number of high-frequency haplotypes. This implies a high turnover of variation and the action of mechanisms that can bring haplotypes to high frequencies relatively rapidly.

A better resolution of the gene genealogy (SIGURGISLASON and ARNASON 2003 ), obtained from a larger fragment of the mtDNA including the TP spacer region (BAKKE et al. 1999 ), increased the variation considerably in terms of both number of haplotypes and nucleotide diversity. The four high-frequency chromosomal types discussed here were also found and overall the structure of a genealogy did not change. Particularly noteworthy is that the distance between the E and the A haplotypes remained one mutation although other branch lengths increased. This is a further indication that haplotype E has arisen recently and its relatively high frequency implies a rapid rise in frequency (SIGURGISLASON and ARNASON 2003 ).

Population differentiation:
Previous studies of mtDNA variation in cod (CARR and MARSHALL 1991A ; PEPIN and CARR 1993 ; CARR et al. 1995 ; ARNASON and PALSSON 1996 ; ARNASON et al. 1998 , ARNASON et al. 2000 ; SIGURGISLASON and ARNASON 2003 ) have not found evidence for significant spatial structuring among populations within a country but weak temporal effects are detected (ARNASON et al. 1998 , ARNASON et al. 2000 ). This study shows clear differentiation on a larger scale across the Atlantic. The AMOVA shows that haplotype frequency differences among countries are a more important source contributing to the differentiation than are mutational differences among haplotypes. Pairwise comparisons show a step in differentiation and gene flow going from Newfoundland or the Baltic to the central countries. However, there is sufficient gene flow to prevent differentiation by random drift alone under an island model. Newfoundland and the Baltic are characterized by haplotypes that are genetically closest together (SIGURGISLASON and ARNASON 2003 ) and nested clade analysis shows that these haplotypes and clades that they represent are overdispersed and geographically displaced relative to other members of the same clade. By inference, population differentiation was due to contemporary forces of isolation by distance but with evidence or long-distance dispersal (and see GULLAND and WILLIAMSSON 1962 ). This indicates a relatively rapid rise in the frequency of these types. Overall, therefore, the differentiation is due primarily to the clines and the high frequencies of the high-frequency haplotypes, in particular the high frequencies of the A and E haplotypes characterizing populations at the opposite ends of the distribution studied here.

The detection of small and statistically significant microspatial structure revealed by microsatellite variation (average F_ST 0.0022–0.0067; RUZZANTE et al. 2001 ; BEACHAM et al. 2002 ; KNUTSEN et al. 2003 ) is not surprising given the mtDNA results, which show that forces operating are not fully contemporaneous across the Atlantic and may therefore be detectable locally. The timescale is important. The microsatellite differentiation, however, is very small and implies considerable gene flow. Such small F_ST is of unknown biological significance in terms of historical divergence and a detailed study of Pan I variation (POGSON and FEVOLDEN 2003 ) failed to find support for the historical hypothesis of deep structure in cod (MOLLER 1969 ; ARNASON and PALSSON 1996 ). The question of temporal stability of population structure shown by microsatellites (RUZZANTE et al. 2001 ; BEACHAM et al. 2002 ) is debatable. For example, inspection of the results (RUZZANTE et al. 2001 ) shows that the spatial structure as determined by F_ST is significantly different from one time period to the next on a decadal scale. This underlines the importance of conducting further studies on temporal variation.

Effective population size:
The shallow genealogy implies a small effective size. Actual population sizes are very large, giving a small N_e/N_a ratio, perhaps as low as 10^–5 and 10^–6 or less. This is low compared to many other organisms (FRANKHAM 1995 ) but similar to high-fecundity marine molluscs and fish (HEDGECOCK et al. 1992 ; TURNER et al. 2002 ). Results of neutral tests (TAJIMA 1989 ; FU 1997 ) do not give unequivocal indications of population expansion, except possibly in Newfoundland. Coalescent methods make more efficient use of data (FELSENSTEIN 1992 ). The results of Markov chain Monte Carlo estimation with fluctuate (KUHNER et al. 1998 ) under a no-growth model were very similar to the nonparametric skyline estimate. The exponential population growth model does not support population expansion. On the contrary the best fit was obtained for a very small and exponentially declining population. The analysis was done by combining the entire sample of sequences from the various populations that are differentiated primarily by frequencies of common haplotypes. The population differentiation is expected to lead to overestimation of effective population size for both fluctuate and skyline estimation. Similarly the estimated negative growth of the best-fit models may also result from the total sample being a mixture from differentiated populations. Under rapid population expansion gene genealogies are expected to be star shaped, that is, having short internal relative to external branches (SLATKIN and HUDSON 1991 ). Population decline has the opposite effect. In the mixed sample most individual sequences coalesce immediately onto the common haplotypes A, E, D, and G. External branches are then short relative to internal branches and the fluctuate estimator returns a negative growth rate. This does not preclude a recent expansion locally, for example, in Newfoundland as the Tajima test suggests. A comparable coalescent estimation for each locality is not done here.

Therefore, the no-growth model is an adequate descriptor of effective population sizes that are small by this method. There is also fair agreement with skyline estimation and conventional estimate based on pairwise differences. Population fluctuation, such as observed in the past 50 years in the Icelandic stock, only modestly reduces effective size and cannot by itself explain the low N_e/N_a ratio. Microsatellite variation shows that F_IS is not different from zero (KNUTSEN et al. 2003 ) and average F_ST is low (see above). Plugging these values into an island-structured metapopulation model (NUNNEY 1999 , Equation 16) gives a negligible effect. If mortality in cod was nonselective such that there was a Poisson distribution of surviving offspring, effective size would be on the order of actual size. Thus neither contemporary population fluctuation nor metapopulation structure can account for the low N_e/N_a ratio. The excess of rare alleles and high diversity observed in RFLP data (POGSON et al. 1995 ) and microsatellites (RUZZANTE et al. 1999 ; BEACHAM et al. 2002 ; KNUTSEN et al. 2003 ) might seem to indicate high effective population sizes at odds with the above. However, these are loci selected for their variation, and less variable microsatellite loci that are not used for population studies exist. Also a combination of high mutation rate, such as at microsatellites, and high actual population sizes makes it likely that one finds rare alleles in population samples because new mutations are realized as individuals in the population. The mitochondrial data show that rare deleterious variants are segregating in the population, which a comparative analysis suggests have not yet been removed by selection.

Atlantic cod is one of the most fecund vertebrates. MAY 1967 found that annual fecundity F was allometric to length L described by the formula log F = 3.42 log L – 0.30. Thus a 100-cm cod could lay 3.5 million eggs, a 150-cm cod 14 million, and a 200-cm cod 37 million. Cod is long lived and can live up to 30 years. Although old and >200-cm fish are very rare, such "mother fish" do exist and are caught by fishermen (e.g., HAEDRICH and FISCHER 1998 ). There are also indications that the coefficient of the allometric equation may increase with age and also that fecundity increases with repeat spawning (MAY 1967 ). Potential lifetime fecundity of cod can best be described as being enormous.

Genotypic variation could influence traits for viability of both eggs and zygotes as well as later developmental stages to adulthood, a form of viability selection. Fecundity selection could also play a role with fitness dependent on genotypes of pairs. Male size also influences male fecundity (BEKKEVOLD et al. 2002 ). In experimental studies of male reproductive competition in spawning aggregations they found, utilizing microsatellite variation, that although most males participated in spawning large males made a disproportionate contribution to fertilization. A study of sperm quality of virgin and repeat spawning males and associated hatching success showed that although success of virgin males was more variable their mean success was the same as the older males (TRIPPEL and NELSON 1992 ). Thus differences between large and small individuals (BEKKEVOLD et al. 2002 ) are likely due to competition and not due to intrinsic differences in sperm quality. Male cod also are long lived and can get very large and presumably have similar fecundity variation as females. Alternatively or concomitantly there may be dilution and turbulence during spawning, leading to sperm limitation (LEVITAN and PETERSEN 1995 ; YUND 2000 ). Thus there are likely to be a multitude of selective forces for successful fertilization on both male and female traits in cod.

Taken together, all these factors translate into a potential for a high variance in fitness and hence a high variance in offspring number because of selection for various traits important for fitness. Thus low N_e/N_a ratio is expected (HEDGECOCK et al. 1992 ; TURNER et al. 2002 ). Genetic loci that are not directly targets for the fitness variation, such as silent mtDNA variation, would be affected by a higher rate of genetic drift.

Five explanations of clines and shallow genealogies:
The important features of the results are the high frequency of the various haplotypes in the different countries; the number and wide geographic distribution of haplotypes spawned off of other haplotypes and a multifurcating constellation genealogy; the trans-Atlantic clines in frequency; the shallowness of the tree of genealogical relationships; and a high apparent rate of homoplasy implying high mutation, that the two haplotypes that are mutationally closest characterize populations at opposite ends of the geographical distribution and the low N_e/N_a ratio. I consider five possible explanations for the features of the data.

First, selection arising from the action of the cellular machinery may be operating on the mutant sites themselves. Sites defining high-frequency haplotypes are synonymous and therefore selection arising from the action of proteins is absent. Natural selection on synonymous sites can, however, arise from the phenotype of the DNA in replication and transcription and from the phenotype of the RNA in translation as the intra- vs. interspecific differences of purine and pyrimidine transitions show. This, however, is an unlikely explanation of the polymorphism for several reasons. First, this is weak purifying selection and not positive selection needed to explain the high frequency of several haplotypes. Also it presumes that by selecting at random a 250-bp fragment of a 16-kb chromosome one finds several selected sites and even balanced polymorphisms due to selection by the cellular machinery.

Second, historical isolation and random drift in refugia might explain the high frequency of different haplotypes in various locations. Assume, for example, that during the Pleistocene the ice sheet created small isolated pockets of suitable habitat for cod and populations of small effective sizes became isolated in various geographic regions. Under these circumstances neutral variation will locally drift to a high frequency or even fixation. With the end of the Pleistocene population sizes expanded and a continuous distribution was established throughout the Atlantic, resulting in a cline because of extensive gene flow among the populations. The cline would, under this scenario, be transient and eventually gene flow would homogenize the frequencies of neutral haplotypes throughout the species range.

Third, this may be a balanced polymorphism. If there are amino acid sites elsewhere on the mitochondrial chromosome that are adapted to local environmental conditions in the waters of each country they could, because of complete linkage of the mitochondrial genome, carry the silent and neutral sites with them to a high frequency in each country. Each of the silent polymorphisms would then simply be a marker for a certain locally adapted chromosome. Thus, under this explanation, haplotype A chromosomes would be adapted to Newfoundland waters, haplotype D chromosomes to Icelandic waters, haplotype E to Baltic waters, and so on. Extensive gene flow establishes clines but local selection is strong enough that it is not swamped by the gene flow. Under this scenario the polymorphism is balanced globally through the interaction of local selection and global gene flow.

Under both the historical and balance explanations deep branches of divergent haplotypes are expected, comparable, for example, to the divergent A and B alleles at the Pan I locus in cod that differ by 19 nucleotides, including radical amino acid changes (POGSON 2001 ). The shallowness of the genealogy is evidence against these explanations. Of all the localities considered the Baltic is the most isolated hydrographically and being the largest brackish water body in the world with occasional influx of sea water from the North Sea replenishing eutrophic conditions, the Baltic also provides a special selective environment for cod (see NISSLING and WESTIN 1991 ; BAGGE and THUROW 1994 ; NISSLING 1994 ). Thus under either explanation the Baltic population would be expected to be most divergent genetically. Yet the Baltic cod is characterized by a high frequency of haplotype E that is mutationally closest to haplotype A (and see SIGURGISLASON and ARNASON 2003 , for resolution of genealogy) that characterizes Newfoundland. On the basis of the totality of the evidence these explanations seem unlikely.

Fourth, there might be frequent selective sweeps of mitochondrial variation, which through linkage have brought haplotypes to high frequencies. Under this explanation mutations that substantially increase the fitness of their carriers frequently arise somewhere on the mitochondrial chromosome. These mutations are adaptive globally and sweep to fixation first locally and then in every population throughout the species range. Thus, perhaps, the A haplotype has recently arisen in Newfoundland and has already almost swept through the Newfoundland population and is spreading eastward, or the E haplotype is sweeping through the Baltic population and is starting to spread to other countries by gene flow and natural selection. Under this scenario a particular chromosome will eventually sweep through the whole population, removing all variation and resetting the stage. The observed polymorphisms would then be the transient phase of current multiple selective sweeps. This explanation can account for the data but the main difficulty is to explain why there would be so much adaptive evolution going on for mitochondrial activity in cod.

Fifth, there may be demic selection of mitochondria. A potential for very high fitness variance and sweepstakes reproduction exists in cod (WILLIAMS 1975 ; HEDGECOCK 1994 ; BEKKEVOLD et al. 2002 ). A mitochondrial type carried by a female winning the reproductive sweepstakes could rapidly increase in frequency. Sweepstakes can be local and in the next generation some offspring of that female again win the local sweepstakes and therefore the frequency of a particular mitochondrial type increases further locally. Some offspring migrate and try their luck elsewhere and so a cline is established via gene flow. The mitochondrial types themselves have nothing to do with the increase in frequency. Which haplotype increases in which locality is haphazard and the memory of the system and hence timescale is a function of interaction between gene flow and variance in net fitness. There is a global polymorphism that depends for stability on an escape in time and space. This assumes high fitness of nuclear genotypes for important fitness traits that are perhaps undergoing selective sweeps and a rapid enough demic spread that the mitochondrial types are not disassociated from the nuclear genotypes. It is possible that some cytoplasmic nuclear interactions would be needed to be invoked to explain the patterns. Alternatively, without any such interactions, a mitochondrial type could get to a high frequency by chance carried by highly fit females. If success in sweepstakes reproduction depends on assembling a highly fit nuclear genotype (WILLIAMS 1975 ) it will most likely carry one of the high-frequency mitochondrial types that would thus increase and decrease in frequency temporally, as seen within countries (ARNASON and PALSSON 1996 ; ARNASON et al. 1998 , ARNASON et al. 2000 ). This would seem therefore to be the most plausible of the five explanations. It is a hypothesis predicting temporal changes. Studies of temporal variation are called for to test it and better resolve the differences between historical and contemporary factors influencing variance in offspring number and effective population sizes.

	ACKNOWLEDGMENTS

I thank the editor and two anonymous reviewers for very helpful and critical comments on the manuscript. This work was supported by grants from the University of Iceland Research Fund and the Icelandic Science Fund.

Manuscript received September 16, 2002; Accepted for publication January 2, 2004.

	LITERATURE CITED

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

ANONYMOUS, 2003 State of marine stocks in Icelandic waters 2002/2003. Prospects for the quota year 2003/2004. Technical Report, Marine Research Institute, Reykjavík, Iceland.

ÁRNASON, E., and J. S. F. BARKER, 2000 Analysis of selection in laboratory and field populations, pp. 182–203 in Evolutionary Genetics: From Molecules to Morphology, edited by R. SINGH and C. KRIMBAS. Columbia University Press, New York.

ÁRNASON, E. and S. PÁLSSON, 1996  Mitochondrial cytochrome b DNA sequence variation of Atlantic cod, Gadus morhua, from Norway. Mol. Ecol. 5:715-724.

ÁRNASON, E., S. PÁLSSON, and P. H. PETERSEN, 1998  Mitochondrial cytochrome b DNA sequence variation of Atlantic cod, Gadus morhua, from the Baltic and the White Seas. Hereditas 129:37-43.[CrossRef][Medline]

ÁRNASON, E., P. H. PETERSEN, K. KRISTINSSON, H. SIGURGÍSLASON, and S. PÁLSSON, 2000  Mitochondrial cytochrome b DNA sequence variation of Atlantic cod from Iceland and Greenland. J. Fish Biol. 56:409-430.[CrossRef]

BAGGE, O. and F. THUROW, 1994  The Baltic cod stock: fluctuations and possible causes. ICES Mar. Sci. Symp. 198:254-268.

BAKKE, I., G. F. SHIELDS, and S. JOHANSEN, 1999  Sequence characterization of a unique intergenic spacer in mtDNA of Gadiformes. Mar. Biotechnol. 1:411-415.[CrossRef]