|Home | About | Journals | Submit | Contact Us | Français|
The ecological theory of adaptive radiation predicts that the evolution of phenotypic diversity within species is generated by divergent natural selection arising from different environments and competition between species. Genetic connectivity among populations is likely also to have an important role in both the origin and maintenance of adaptive genetic diversity. Our goal was to evaluate the potential roles of genetic connectivity and natural selection in the maintenance of adaptive phenotypic differences among morphs of Arctic charr, Salvelinus alpinus, in Iceland. At a large spatial scale, we tested the predictive power of geographic structure and phenotypic variation for patterns of neutral genetic variation among populations throughout Iceland. At a smaller scale, we evaluated the genetic differentiation between two morphs in Lake Thingvallavatn relative to historically explicit, coalescent-based null models of the evolutionary history of these lineages. At the large spatial scale, populations are highly differentiated, but weakly structured, both geographically and with respect to patterns of phenotypic variation. At the intralacustrine scale, we observe modest genetic differentiation between two morphs, but this level of differentiation is nonetheless consistent with strong reproductive isolation throughout the Holocene. Rather than a result of the homogenizing effect of gene flow in a system at migration-drift equilibrium, the modest level of genetic differentiation could equally be a result of slow neutral divergence by drift in large populations. We conclude that contemporary and recent patterns of restricted gene flow have been highly conducive to the evolution and maintenance of adaptive genetic variation in Icelandic Arctic charr.
Evidence for the important role of natural selection in the generation of diversity, such as the repeated pattern of phenotypic divergence connected to differential exploitation of discrete habitats, is overwhelming (Schluter, 2000). Although such divergence must result from selection acting on the available genetic variation (Barrett and Schluter, 2007), the availability, quantity and other aspects of the nature of this variation are determined by the other evolutionary forces: mutation, and especially at microevolutionary scales, drift and migration. In particular, natural selection is expected, almost universally, to erode genetic variation in populations, and thus the existence of the raw material for evolutionary change depends primarily on non-selective processes. Importantly, the non-selective evolutionary forces can all interact with selection in constraining or facilitating modes (Garant et al., 2007) and can vary both spatially and temporally. Natural selection will vary with ecological conditions; migration will vary with patterns of genetic connectivity; and the consequences of drift and mutation will be highly dependent on variation in local population size (Kimura, 1983; Caballero, 1994). Thus the opportunity for adaptive diversification is necessarily geographically structured by spatial variation in the efficacy and interactions of the different evolutionary forces.
It is particularly important to consider the evolution of local adaptation within species in the context of prevailing patterns of gene flow. Gene flow can both constrain and facilitate adaptive divergence (Garant et al., 2007). The constraining potential of gene flow has a strong and straightforward theoretical basis (Slatkin, 1973; Felsenstein, 1976; Lenormand, 2002), and has been demonstrated by the inverse relationships between gene flow and adaptive divergence (Hendry et al., 2002; Garant et al., 2007). This pattern has been described in several species including threespine stickleback, Gasterosteus aculeatus, in which gene flow from lake to river morphs appears to constrain morphological divergence (Hendry and Taylor, 2004; Moore and Hendry, 2005), and similar patterns in ecologically relevant traits in other taxa have been reported as well (for example, Wood and Foote, 1996; Smith et al., 1997; Lu and Bernatchez, 1999; Fraser and Bernatchez, 2005). However, low or intermediate levels of gene flow can also be beneficial to adaptive divergence (Gomulkiewicz et al., 1999; Doebeli and Dieckmann, 2003; Garant et al., 2007). Gene flow can maintain genetic diversity, thereby increasing adaptive potential (Swindell and Bouzat, 2006), and can counteract inbreeding depression especially in small populations (Ingvarsson and Whitlock, 2000; Ingvarsson, 2001). Thus, our understanding of the processes responsible for adaptive divergence in any species will be incomplete without consideration of the role of gene flow and, most specifically, the genetic connectivity among natural populations.
In addition to spatial variation, the extent of adaptation, diversity, the role and efficacy of selection, and the influence of other evolutionary forces undoubtedly change over time. For example, in threespine stickleback, recent changes in water clarity of Enos Lake have reduced the potential for assortative mating between nascent species that were adapted to different environments within the lake (Taylor et al., 2006). This has caused an increase in gene flow, and the ability of natural selection to maintain the adaptive diversity has not been sufficient to overcome the recent increase in the exchange of alleles between the two types. Along similar lines, Butlin et al. (2008) present a hypothesis for the origin of apparent parallel patterns, in which connectivity and selection both act constructively, but at different times. Under their model, alleles that are advantageous in a particular type of habitat arise by mutation only once, but become geographically widespread because of rare migration events, and are subsequently independently selected and introgressed into the genomes of organisms in widespread locations.
In many species, the distribution of genetic diversity may be much more the product of past processes than current forces. For example, a large diversity of mitochondrial haplotypes are distributed widely among lake trout Salvelinus namaycush throughout the central portion of their range in North America, and this distribution is a result of the mixing of divergent lineages in proglacial lakes that no longer exist (Wilson and Hebert, 1998). Such mixing of diverent lineages would likely have generated diversity, and importantly, this feature of lake trout genetics would not be evident from consideration of current spatial arrangements. Similarly, effects of transient patterns of postglacial connectivity on the distribution of genetic diversity are common at various scales in other fish species (for example, Poissant et al., 2005), and in more mobile organisms such as birds, vicariance during the Pleistocene is often the cause of major patterns of within-species variation (Avise and Walker, 1998).
Details of historical demographic processes (for example, dispersal from origins or refugia, connectivity/gene flow, bottlenecks and so on) are largely unknown for most species, and consequently, inference of the cause of existing patterns of adaptive variation is often difficult. For example, repeated phenotype–environment matching may result from multiple evolutionary origins of locally adapted phenotypes (Svardson, 1979) or may result from the widespread dispersal of apparently locally adapted morphs, after having evolved once or relatively few times (Behnke, 1972). However, if the occurrence of subsequent gene flow cannot be ruled out, it is very difficult to distinguish between these alternate hypotheses (Ferguson and Mason, 1981; Schluter, 1996; Volpe and Ferguson, 1996; Taylor, 1999; Butlin et al., 2008).
The evolution of phenotypically distinct morphs is common in north-temperate fish species. Reports exist from several families, including the Gasterosteidae (threespine stickleback; McPhail, 1993), Osmeridae (rainbow smelts Osmerus mordax; Taylor and Bentzen, 1993), Centrarchidae (pumpkinseed sunfish Lepomis gibbosus; Robinson et al., 1993) and especially the Salmonidae, in which examples are available from a number of species, including lake whitefish Coregonus clupeaformis (Fenderson, 1964), lake trout (Blackie et al., 2003) and brown trout Salmo trutta (Ferguson and Mason, 1981). This has been particularly well documented in Arctic charr, Salvelinus alpinus, in Iceland (Snorrason and Skúlason, 2004). In Icelandic lakes, the co-occurrence of phenotypically distinct morphs of Arctic charr is common (Gíslason et al., 1999; Wilson et al., 2004). This variation clearly correlates with the ecological and physical diversity of habitats, and it has been suggested that the special attributes of the neovolcanic zone, with extensive lava fields and springwater systems together with the paucity of the freshwater fish fauna, may have been a key causal agent of this observed diversity (Skúlason et al., 1999; Snorrason and Skúlason, 2004).
A striking example of intralacustrine diversity in Arctic charr is present in Lake Thingvallavatn (Figure 1), where distinct piscivorous, pelagic (planktivorous) and large and small benthic morphs exist (Snorrason et al., 1989). Natural history provides suggestive evidence for the adaptive nature of this variation, in which the diversity in a range of phenotypic traits, including life histories, behavior and morphology covary closely with the ecological features of each morph's niche (Sandlund et al., 1987; Jonsson et al., 1988; Skúlason et al., 1989a, 1989b; Malmquist et al., 1992). Common-garden experiments have indicated that the phenotypic differences are the result of both genetic differences (Skúlason et al., 1989a, 1993) and genetically based differences in putatively adaptive plastic responses to different environments (Parsons et al., 2010). At a larger spatial scale, morphs that are phenotypically diverged along similar axes inhabit other lakes (Gíslason et al., 1999; Wilson et al. 2004), and a derived small benthic morph associated with and apparently adapted to benthic environments in small springs is particularly widespread (Kristjánsson, 2008).
Although the signature of natural selection is highly evident in the repeated patterns of divergence and repeated phenotype–environment matches in Icelandic Arctic charr populations, including both intralacustrine morphs and the widespread small benthic morph, we know very little about the extent to which different evolutionary forces have contributed to the development of this variation. For example, we do not know whether or not the different morphs evolved multiple times in situ, or whether their widespread distribution is attributable more to dispersal of the different morphs. Regardless of the process that initially created the widespread, apparently repeated locally adapted morphs, selection has undoubtedly had some role in the maintenance of this diversity. However, we do not know to what extent this inferred selection has been necessary to counteract homogenizing effects of gene flow, or whether the populations have been relatively isolated for an extended period of time.
Our primary goal herein was to determine the extent of connectivity among populations of Arctic charr throughout Iceland and within Lake Thingvallavatn to determine the relative importance of gene flow and selection in the maintenance of adaptive spatial variation in phenotype. Second, we sought to determine the extent to which variation in connectivity throughout the Holocene has shaped current patterns of neutral variation to make inferences about the relative extents to which gene flow and selection have contributed to generation of contemporary patterns of neutral molecular and adaptive phenotypic variation. In Lake Thingvallavatn, we compared observed levels of genetic differentiation with distributions of population-genetic statistics generated under stochastic models to determine what historical patterns of connectivity might plausibly have generated currently observed distributions of genetic variation. At the Iceland-wide scale, we used waterway and overland geographic distance in conjunction with the occurrence of the small benthic morph to predict variation in genetic differentiation among populations. This allowed us to infer the extent to which patterns of genetic variation are consequences of processes during postglacial colonization, relative to current patterns of connectivity. Consequently, we were able to determine whether or not the small benthic morph has had the opportunity to evolve multiple times, or whether its widespread distribution is more likely to be the result of dispersal after it first evolved.
Icelandic Arctic charr populations vary extensively in morphology, behavior and life history (Snorrason and Skúlason, 2004; Kristjánsson, 2008). As mentioned previously, Lake Thingvallavatn contains four morphs, and in the current study, we focus on the small benthic and pelagic morphs. The small benthic morph is dark in color, has relatively few gill rakers, stocky body, long pectoral fins, blunt snout and short lower jaw. These fishes live mainly in the stony littoral zone and feed on benthic invertebrates, especially on the gastropod Lymnaea peregrea. In contrast, planktivorous charr are silvery, have higher number of gill rakers, a fusiform body, smaller pectoral fins, a pointed snout and longer lower jaw. They feed in the water column on zooplankton and ascending chironomid pupae (Sandlund et al., 1987; Snorrason et al., 1989; Malmquist et al., 1992). The morphs differ extensively in their life history characteristics. Small benthivorous charr mature at 2 (males) to 4 years (females) and at a minimum fork length of 7.2cm (males), and planktivorous charr mature at 3–5 years and at the minimum fork length of 15.2cm (Jonsson et al., 1988). The small benthic and pelagic morphs both spawn in the stony littoral zone around the lake at depths from 0 to 8m and overlap in spawning time (September–October; Skúlason et al., 1989b).
Lake Thingvallavatn is Iceland's largest lake with surface area of 83.7km2. It is located in the axial rift zone of southwestern Iceland at 100.7m above sea level. As with most of the watershed, the area north of Lake Thingvallavatn is covered by postglacial (formed after 10000ybp) lavas, and the major part of the water entering the lake is spring water, pouring through rifts on the north and east shores (Adalsteinsson et al., 1992). A lake about the present lake's size was formed after the last glaciation ~10000 years ago (Sæmundsson, 1992). Lava flows and isostatic rebound, together with formation of waterfalls on the outlet river, probably caused isolation of Lake Thingvallavatn soon thereafter, thereby cutting off possibilities of immigration from downstream populations of Arctic charr. High tectonic activity with rift formations and eruptions may have prompted transient isolations of small water bodies (for example, in rifts) during this period, thus potentially temporarily closing off small populations of Arctic charr. The time period between 10000 and 12500ybp, before and during the last transient advance of ice, is of considerable interest. Early in this period, sea levels were very high in southwestern Iceland (Norðdahl and Pétursson, 2005), and for a brief geological period, the sea may even have flooded the Lake Thingvallavatn depression. The extent of maximum ice cover during the last transient glacial advance (between 10000 and 11000ybp) is not well known for the Lake Thingvallavatn area. There is, however, strong evidence that small glacial lakes existed between the glacier and the Hengill mountain (that is, within the present Lake Thingvallavatn depression; Sæmundsson, 1992). If these lakes were remnants of the lake/fjord that was formed earlier, the possibility exists that Arctic charr may have persisted there throughout the last transient glacial advance. Thus, Lake Thingvallavatn has likely contained Arctic charr for the last 10000 years. However, the exact time of the first colonization is unclear, and whether or not multiple colonizations or (micro) allopatric events have contributed to the contemporary genetic structure of the lake's Arctic charr populations is unknown.
The widespread occurrence of a small benthic morph similar to the small benthivorous charr of Lake Thingvallavatn has recently been documented in the neovolcanic zone of Iceland (Sigursteinsdóttir and Kristjánsson, 2005; Kristjánsson, 2008). These fishes typically inhabit spring habitats within or near lava fields, which in some cases appear to be physically isolated from other water bodies. These small fishes tend to be similar in terms of life history and morphology. Typically, they exhibit small size at maturity, low growth rate after maturity and a juvenile appearance with blunt snout and cryptic coloration, for example they retain their parr marks (Kristjánsson, 2008). Broadly, this phenotype is consistent with adaptation to small, spring water habitats, often with low temperature the year round and a spatial structure that favors small size, and with high fish density (Snorrason et al., 1994; Snorrason and Skúlason, 2004; Kristjánsson, 2008). Despite subtle differences in diet and morphology among populations of this morph, they all differ markedly in these characters from typical phenotypes of Arctic charr (Sigursteinsdóttir and Kristjánsson, 2005; Kristjánsson, 2008).
We collected small benthic charr from 20 locations in several drainages in the northeast and southwest of Iceland (Table 1 and Figure 1). Additionally, we sampled eight reference populations. We chose reference populations on the basis of geographical proximity to areas from which we had sampled small benthic populations. Our reference populations of charr are characterised by larger adult body size, and were obtained primarily from anadromous populations. When proximate anadromous populations did not exist, we obtained reference samples from large-bodied, morphologically unspecialized lake-resident populations. Consistent sampling of a single alternate morph or of a range of morphs from each geographic area was not possible because of unevenness of the spatial distribution of the different morphs. We captured all fish by electrofishing or by gill netting between 2003 and 2006. Sample size at each location varied between 38 and 46 individuals. We sacrificed each fish, removed the right pectoral fin and stored it in 96% ethanol for subsequent molecular analyses.
We captured 503 small benthic and pelagic charr in October 2005 from five locations around Lake Thingvallavatn (Figure 1, inset) by gill netting. The two morphs were classified according to Snorrason et al. (1994), see above. All fishes were mature with ripe gonads. We also obtained a sample of 50 large benthic charr from Ólafsdráttur (Figure 1), their principal spawning ground within Lake Thingvallavatn, for comparison. We obtained and stored samples of fin tissue for each of these fishes as described above. The large benthic reference sample and the sample of small benthic charr from Ólafsdráttur were included in the Iceland-wide sample set, described above.
We obtained genomic DNA from all fin tissue samples using the GenElute Mammalian Genomic DNA Miniprep Kit (Sigma Corp, St Louis, MO, USA). We genotyped all samples at nine microsatellite loci (bracketed values are Genebank accession numbers, when not citations): BHMS206 (AF256680), BHMS417 (AF256752), Bx079862 (BX079862), Omi127 (AB105850), Omi187 (AB105857), OMM1236 (AF470016), One11ASC (Scribner et al., 1996), SalD25SFU (McGowan et al., 2004) and Sco19SFU (Taylor et al., 2001). Omi187 amplifies two loci in Icelandic Arctic charr (a common occurrence in salmonids due to their ancestral autotetraploidy; Allendorf and Thorgaard, 1984), both of which are polymorphic in samples from Lake Thingvallavatn and in linkage equilibrium with one another (see below). Consequently, our marker set consists of 9 loci for analyses conducted on the Iceland-wide scale and 10 loci for the within-Lake Thingvallavatn analyses. We amplified each locus separately by PCR following the protocol described in Wilson et al. (2004) and using predetermined locus-specific annealing temperatures (Appendix Table A1). One primer for each locus was fluorescently labeled to allow subsequent visualization using an Hitachi FMBIOII fluorescence imaging system (Hitachi, San Francisco, CA, USA) following electrophoresis in 6% polyacrylamide gels. All alleles were sized relative to multiple 350-TAMRA lane standards (Applied Biosystems, Carlsbad, CA, USA) that we ran in each gel.
We conducted a range of analyses to characterize our loci and to confirm the reliability of our marker set for subsequent analyses, both at the intralacustrine and at the Iceland-wide scales. First, we calculated indices of genetic diversity, including the number of observed alleles, expected heterozygosity, and allelic richness using MSA 4.05 (Dieringer and Shlötterer, 2003). We tested for deviations from Hardy–Weinberg proportions and linkage disequilibrium using the permutation tests implemented in GENETIX (Belkhir et al., 2003). To accommodate multiple tests, we applied Bonferroni corrections (Rice, 1989) at the level of populations for tests of observed and expected heterozygosity, and at the level of pairs of loci for test of linkage disequilibrium.
We conducted two complementary analyses to determine the extent to which genetic differentation among charr populations is related to morph differences, both in the context of spatial variation and relative to the influence of genetic connectivity. We calculated hierarchical F-statistics (Weir, 1990; Yang, 1998) to assess the potential importance of morph classification, relative to spatial variation. We used morph (small benthic versus reference), region (northeast versus southwest) based on the two large-scale spatial clusters of our samples, watershed and population as nested hierarchical levels. We constructed three hierarchies based on these levels, with all possible combinations of region, watershed and morph, and always with population at the lowest level. For each hierarchy, we calculated the variance (F) associated with each level using the function varcomp.glob in the R package hierfstat (Goudet, 2005). We calculated the statistical significance of variance explained at each level using the permutation algorithms implemented in the functions test.between (highest levels), test.between.within (intermediate levels) and test.within (lowest, that is, population, level), implemented in the same package. To visualize similarity among populations, we generated a phenogram based on Cavalli–Sforza's genetic distance and the neighbor-joining algorithm implemented in the package PHYLIP (Felsenstein, 1989).
To treat spatial variation as a continuous variable in our analyses of the relationship between morph and genetic subdivision, we adopted a multiple regression-like Mantel test approach similar to that described by Smouse et al. (1986). To separate geographic effects acting on different scales, we modeled a predictor of Cavalli–Sforza's (Cavalli-Sforza and Edwards, 1967) pairwise genetic distance E(Dij) as
where i and j are index populations, δr is an indicator variable with a value of 1 when populations i and j are in different regions and 0 otherwise, δm is an indicator variable with a value of 1 when populations are of different morphs and 0 otherwise. A is a coefficient predicting genetic distance between regions, B is a coefficient relating expected genetic distance to Gij, the geographic distances between populations, and C is a coefficient relating morph differences to genetic distance. Equation (1a) thus models morph as a predictor of genetic distance at the scale of the entire analysis, and equation (1b) models morph as a predictor of genetic distance only within region. In equations (1a and 1b), B is only a predictor of differentiation within regions; as a consequence of the fact that the middle terms on the left hand side are zero when δr=1. Consequently, large-scale geographic variation in allele frequencies will not drive inferences of isolation by distance at smaller scales. We calculated the correlation between matrices of observed genetic distance and E(Dij) values using the function mantel.rtest in the R package ade4 (Dray and Dufour, 2007), and obtained the values of the coefficients A, B, and C that maximize the correlation between the matrices using the modified Newton–Raphson algorithm implemented by the function nlm in the R package base (R Development Core Team, 2005). We tested both Euclidean overland distance and waterway distance for values of Gij, whereby testing both direct distance and waterway distances as predictors of genetic distance.
We tested for and characterized genetic subdivision within and among the small benthic and pelagic morphs of Arctic charr in Lake Thingvallavatn. We calculated Weir and Cockerham's (1984) estimator of Fst, θ, using FSTAT (Goudet, 1995) for (a) global differentiation including the large benthic reference sample, (b) global differentiation including only the small benthic and pelagic samples and (c) all pairwise comparisons of samples. We obtained the statistical significance of all θ estimates from 1000 bootstrap samples. Next, we conducted analyses of molecular variance to determine more quantitatively the extent to which genetic variance is associated with morphs, populations and individuals within populations. We implemented the analyses of molecular variance with the software ARLEQUIN (Excoffier et al., 2005). Finally, we generated a neighbor-joining phenogram for samples collected within Lake Thingvallavatn, as described above.
To explore which historical scenarios could be responsible for the current distribution of microsatellite allele diversity within and among morphs of Arctic charr in Lake Thingvallavatn, we developed two contrasting plausible models of the evolution of allele frequency differences between the small benthic and pelagic morphs of Arctic charr in the lake. The first model (Figure 2a) represents a scenario of sympatric divergence, whereas the second model involves sympatry throughout most of the Holocene, but preceeded by a transient allopatric phase (Figure 2b) associated roughly with the early postglacial period during which sea level fluctuations and drainage rearrangements following the last transient glacial advance were occurring in Iceland. The first model is thus intended as a simple scenario, where the present Arctic charr fauna of Lake Thingvallavatn are derived from a single invasion, and therefore the differentiation between small benthic and pelagic morphs would have to have happened entirely in sympatry. The time frame, in generations, is determined roughly by the time since deglaciation of the Lake Thingvallavatn area (Sæmundsson, 1992; and see above) and the generation times of Arctic charr (Jonsson et al., 1988). The second model is intended to be a contrasting scenario, wherein the complex geological history of the area may have provided opportunities for the initiation of genetic and phenotypic differentiation to occur in allopatry. This could result from either differentiation in small isolated waterbodies or multiple invasions or some combination of such processes. The allopatric period in the second model is intended to correspond to the period between about 12000ybp, when a lake first formed in the Lake Thingvallavatn area, and when the lake that has since maintained approximately its current size and has been isolated from the sea formed about 10000ybp (see ‘Study system' section, above).
For large ranges of population size and migration rates between morphs, we conducted coalescent simulations based on these two historical scenarios. Briefly, coalescent simulation is a flexible technique whereby (typically large numbers of) hypothetical ancestries of a contemporary sample are generated according to the Fisher–Wright model of an idealized population, given demographic parameters specifying deme number, local effective population sizes, migration rates, as well as potential variation in these parameters, depending on historical hypotheses (Rosenberg and Nordborg, 2002). Distributions of expected genetic variation, given parameters specifying the mutation model and mutation rates, can then be generated, given the sample of hypothetical geneologies. Observed patterns of genetic variation within and among populations can then be compared with models of historical demographic parameters to determine what scenarios can plausibly have generated the parameters or statistics obtained from empirical data.
More specifically, we simulated two contemporary demes with population sizes varying between 103.5 and 106.5. These population sizes bracket the range of estimated contemporary population size (Snorrason et al., 1992) and also include a much lower long-term effective population size than currently exists. We simulated a wide range of migration rates, including values of individual migration rate, m between 10−3 and 10−8. In population-wide terms, this corresponds to a very wide range of total numbers of migrants, that is, Nem from very near zero to many migrants between morphs per generation. We simulated these parameters back in time either 4000 (sympatric scenario; Figure 2a) or 3500 (micro-allopatric scenario; Figure 2b) generations. In the sympatric scenario, we then simulated colonization by movement of all existing lineages into a single deme of size N=1000. In the micro-allopatric scenario, we simulated a reduction of population size in each deme to N=1000 for 500 generations, without migration between morphs, followed by movement into a single deme of size N=1000. Each coalescent simulation generated hypothetical geneologies for 10 loci, with sample sizes equal to the number of samples we collected from each morph in Lake Thingvallavatn. We simulated microsatellite genotypes on each simulated geneology. Because we do not know the mutation rates of the loci we genotyped in this study, we sampled hypothetical mutation rates from the distribution 10δ(−3.8, 0.3), where δ represents a normal distribution with the specified mean and standard deviation. We simulated a stepwise mutation model. This generated a distribution of mutation rates similar to the ranges of estimated microsatellite mutation rates in fishes that have been reported in the literature (especially in the zebrafish map, in which the estimate is based on many loci; Shimoda et al., 1999). We sampled this distribution of mutation rates independently for the generation of microsatellite genotypes for each locus in each coalescent simulation. We conducted 100 replicates of each of the 143 combinations of migration rate and population size for each of the historical scenarios. We conducted all coalescent simulations with simcoal2 (Laval and Excoffier, 2004).
We calculated θ and He for the simulated microsatellite data for each simulation scenario. We then calculated the means of θ and of He for each scenario and the proportions of simulations within each scenario that generated values of θ and of He with absolute deviations from the mean that were greater than the deviation of the observed data from the means. Essentially, we quantified how well the observed data matched the distribution of simulated data for each scenario. To obtain a simple overall metric of the fit between the simulated and observed data for each of the 286 simulated scenarios, we multiplied the numerically obtained probabilities of obtaining the observed data for each of the two statistics.
Within the populations that we sampled from across Iceland, the mean number of observed alleles per locus ranged from 3.67 to 11.44 and observed heterozygosity ranged from 0.377 to 0.696 (Table 1). Within samples from Lake Thingvallavatn, the mean number of alleles per locus ranged from 7.90 to 8.90 and observed heterozygosity ranged from 0.530 to 0.644. Following Bonferroni correction, we detected significant deviations of observed from expected levels of heterozygosity in two populations (S-A-04, N-A-25) and both were associated with a deficit of heterozygotes. This is a common finding in Arctic charr populations, and is generally attributed to cryptic population structure within samples (for example, Wilson et al., 2004). Pooling all within-population tests among loci for linkage disequilibrium, no significant associations of alleles between loci were detected. We detected neither significant differences between observed and expected heterozygosities nor significant linkage disequilibria in samples collected from Lake Thingvallavatn.
All populations of Arctic charr that we included in the Iceland-wide analyses are significantly differentiated from one another, with an overall FST estimate (θ) of 0.245 (95% confidence interval (CI): 0.182–0.311). Mean pairwise FST estimates among small benthic and reference samples are very similar, with values of 0.242 and 0.244, respectively (Appendix Table A2). Despite this high level of allele frequency variation among populations, patterns of similarity among populations are not well resolved based on the application of the neighbour-joining algorithm to estimated genetic distance (Figure 3a), in which many populations are included in a large central node in the unrooted phenogram. This main unresolved group includes both small benthic populations and reference populations, as well as samples from both the northwestern and southwestern regions. Resolved structure in the phenogram largely reflects geographical variation, in which all nodes with >50% bootstrap support contain only populations from the same region. The one case in which multiple samples of the same morph group together is a case in which multiple small benthic populations exist in close proximity (S-SB-11 to S-SB-15), without a nearby reference sample.
Waterway distance is not a significant predictor of genetic distance among populations of Arctic charr within the northeastern and southwestern regions of Iceland. Direct overland distance, however, is a highly significant predictor of genetic distance, and substantially improves the predictive power of the models described by equations (1a and 1b) (Table 2). Specifically, modeling only an effect of region generates a correlation (r2) of 0.127 with estimated genetic distance among samples. Inclusion of waterway distance within regions improves the model fit only to r2=0.149 (P=0.09, one-tailed bootstrap test). Inclusion of direct overland distance, in addition to region, improves the model fit from r2=0.127 to r2=0.270 (P<0.001, one-tailed bootstrap test). We therefore conducted subsequent tests for effects of morph using only overland distance as a geographic predictor.
Inclusion of morph as a predictor of genetic distance produced only small increases in the model fit (Table 2). When we fit morph dissimilarity as a predictor of genetic distance at both within- and between-region levels, the model improvement was small and non-significant. When we fit dissimilarity as a predictor of genetic distance at only the within-region level, we found a significant but small negative relationship, that is, samples from different morphs were slightly, but highly significantly, more genetically similar than samples of the same morph, conditional on geographic distance (P>0.999 for the one-tailed bootstrap test, in which the alternate hypothesis is that morph dissimilarity is positively associated with genetic distance; Table 2).
The hierarchical analyses based on F-statistics corroborate the general finding that morph is a poor predictor of genetic distance. Variances associated with morph were small and slightly negative, regardless of the hierarchial level at which morph was included (Table 3a). Conditional on watershed, region was a minor predictor of the distribution of genetic variation as well (Table 3a), with most of the genetic variation being attributable to the levels of watershed and population.
The overall estimate of FST (that is, θ, Weir and Cockerham, 1984) treating all samples from Lake Thingvallavatn as separate populations is 0.039 (95% CI: 0.012–0.059). Pairwise FST estimates are generally an order of magnitude lower than we observed among the samples collected throughout Iceland (Appendix Table A3). Pairwise θ values between the large benthic morph and the two focal morphs are larger than those among or between samples of the small benthic and pelagic morphs (Appendix Table A3). Pooling samples within morphs, that is, one population for each morph, and excluding the large benthic reference morph samples, θ is 0.036 (95% CI: 0.008–0.067). Within morphs, θ values are 0.003 (95% CI: 0.001–0.005) and 0.006 (95% CI: 0.003–0.012) in the small benthic and pelagic morphs, respectively. In the hierarchical analysis of molecular variance, a small (3.9% Table 3b) but significant (P<0.001) proportion of the total genetic variance is attributable to the morph level (Table 3), and the variance among locations is only statistically significant when nested within morphs (Table 3bii). This finding corroborates the structure of the neighbour-joining phenogram (Figure 3b), in which both morphs cluster separately.
The ranges of hypothetical historical scenarios for which we conducted coalescent simulations generated ranges of θ that spanned the observed values (Figures 4a and b) and ranges of He that were generally close to the observed values, except at the smallest simulated population sizes (Figures 4c and d). In general, population sizes smaller than 104 resulted in the evolution of more genetic differentiation and less genetic diversity than is observed in small benthic and pelagic Arctic charr in Lake Thingvallavatn. In the sympatric model (Figures 2a and 4a, c, e), the range of simulation scenarios with population sizes larger than about 105 generated substantially less genetic differentiation between morphs than we observed. Conversely, the entire range of large population sizes, except at the highest contemporary migration rates, generated approximately the observed level of genetic differentiation in the micro-allopatric scenario (Figures 2b and 4b, d, f). In both the sympatric and micro-allopatric scenarios, some combinations of relatively small popualtion sizes and relatively high migration rates generated similar θ and He values to those observed in the real data of 0.036 and 0.60, respectively. However, these scenarios generated similar patterns in substantially smaller fractions of simulations than scenarios with larger population sizes and smaller contemporary migration rates.
Recent and contemporary gene flow among populations of Arctic charr throughout Iceland has been and is generally very restricted. Across Iceland, we found substantial genetic subdivision, indicating that connectivity is low. Furthermore, direct overland distance predicts genetic distance far better than does waterway distance (Table 3), suggesting that the genetic signatures of processes during the colonization of Iceland by Arctic charr early in the Holocene have not been masked by subsequent gene flow. Additionally, the general lack of resolution of the distance-based phenogram (Figure 3a) despite substantial genetic subdivision (Table 3; see also Appendix Table A2) is consistent with substantial, or in some cases total isolation following the initial postglacial colonization of Iceland by Arctic charr. Given this low genetic connectivity, it is unlikely that the discordance between patterns of phenotypic and genetic variation (Tables 2 and and3)3) results from recent gene flow masking a common origin of the small benthic morph. In other words, had the small benthic morph evolved once and then dispersed widely, the signature of genetic similarity among phenotypically similar forms would remain. Thus, the extreme hypothesis of a single origin of the small benthic morph with subsequent widespread dispersal seems untenable. It is much more probable that the small benthic morph has evolved multiple times throughout the range that we have sampled.
Generally, it seems likely that much of Iceland was colonized by Arctic charr immediately following deglaciation, when the isostatic depression of the region caused relatively high sea levels (Norðdahl and Pétursson, 2005). With rapid isostatic rebound, many Icelandic Arctic charr populations would have become and remained highly isolated. This is consistent with our observation of a large central group in the neighbour-joining phenogram based on genetic distance (Figure 3), despite high genetic subdivision, and is also consistent with the lack of isolation by current waterway distance (Table 2).
It is tempting to infer that adaptive phenotypes have arisen multiple times when particular phenotypes are repeatedly associated with common environments and when geographically proximate but phenotypically dissimilar populations are most genetically similar (Skúlason et al., 1996; Volpe and Ferguson, 1996). However, this interpretation must always be approached with caution because alternative hypotheses exist. The same geographic patterns in adaptive phenotypic variation and genetic variation are compatible with a single evolutionary origin of different morphs, followed by widespread dispersal, but with subsequent gene flow. Such gene flow could cause the observed genetic patterns at marker loci, and indeed, most authors have been very careful to account for this or similar scenarios (Schluter, 1996; Volpe and Ferguson, 1996). In fact, intermediate scenarios should be acknowledged, and therein the truth probably lies. In our study of the contemporary and recent demographics of the small benthic morph of Arctic charr in Iceland, we have been able to largely rule out the latter extreme hypothesis. This illustrates the importance of understanding the demographic context in which adaptive processes occur. By doing this, we have not only been better able to understand the role of selection relative to the other evolutionary forces, but the detailed consideration of patterns of connectivity also provided the means to distinguish between the two alternate hypotheses about the nature of the adaptive origin of the small benthic morph.
Although the lack of an overall pattern of structuring by current waterway arrangements combined with the lack of structuring by morph strongly supports a scenario of multiple independent origins, it does not necessarily follow that the apparently adaptive phenotypes are independently derived in every small benthic population. For example, in the River Ölfusá drainage in the south, there is some local clustering by morph, where the reference populations from Alftavatn and Grimsnes (S-A-01 and S-A-04; Table 1, Figure 1) cluster together (Figure 3a), but both are anadromous populations entering the sea by the River Ölfusá. Similarly, some of the small benthic populations are geographically clustered in a small area, where no proximate sample of a different morph was available for comparison (especially samples S-SB-11 through S-SB-15; Figures 1 and and3a).3a). Hence, although inferences are not possible regarding connectivity among morphs in all areas of Iceland, overall patterns of connectivity have generally been such that adaptive differentiation has not been constrained by the homogenizing effects of gene flow.
We cannot rule out, however, the possibility that some low level of gene flow has contributed to the adaptive phenotypic evolution of the small benthic morph. For example, it is probably naive to conclude that gene flow did not contribute to the origin of the small benthic morph, simply on the basis that this gene flow must not have been recent and must not have occurred along the current arrangement of waterways. During the late Pleistocene and early Holocene, the drainage arrangements of Iceland, just like in other recently deglaciated landscapes (for example, Fulton, 1989), were highly dynamic (Norðdahl and Pétursson, 2005). We do not know enough about the specific nature of drainage rearrangements to form specific hypotheses about the genetic consequences of geomorphological processes during this period. Hence, although movement of alleles that confer adaptation to benthic environments could have happened during this period, promoting the subsequent evolvability of various morphs, it seems likely that the bulk of the evolution to the locally adapted small benthic state has happened in situ.
In the vicinity of Lake Thingvallavatn, we do know enough about the postglacial geography (Sæmundsson, 1992) to construct useful, if simplistic, historical hypotheses about the evolution of phenotypic diversity and adaptation in Arctic charr. In and of itself, the relatively low level of genetic differentiation among the small benthic and pelagic morphs in Lake Thingvallavatn could be taken as evidence for gene flow. In this case, one would infer substantial potential for non-selective forces to constrain adaptive differentiation, especially relative to the general pattern of isolation that we found at the Iceland-wide scale. This is a common interpretation of low levels of genetic differentiation, but depends on the metapopulation being at mutation-drift equilibrium. However, populations of small benthic and pelagic Arctic charr in Lake Thingvallavatn are most likely very large in comparison with most other populations in Iceland, and certainly many orders of magnitude larger than most populations of the small benthic morph. Importantly, the temporal persistence of non-equilibrium patterns of neutral genetic variation is positively related to population size (Whitlock, 1992). In other words, in the absence of gene flow, genetic differentiation will evolve very slowly, if population sizes are large.
The modest level of genetic subdivision between small benthic and pelagic charr in Lake Thingvallavatn is in fact consistent with very low levels of gene flow throughout most of the Holocene. This conclusion is conditional on relatively large population sizes throughout this period. Given the various simplifying assumptions in our historical models, the population sizes of each morph required to generate the observed levels of subdivision and diversity are on the order of 104–105 individuals in the purely sympatric scenario and over 104 individuals in the micro-allopatric scenario (Figure 4). These are in fact modest population sizes for such a large lake, and are indeed very much in the lower end of the range of estimated contemporary population sizes (Snorrason et al., 1992). Long-term variation in population sizes can substantially decrease genetically effective population sizes, because long-term effective population size is proportional to the harmonic mean of population sizes over time (Caballero, 1994). Consequently, we can neither exclude the possibilities that long-run effective population sizes are large nor small. However, despite the ongoing geological instability of Iceland, including activity in and around Lake Thingvallavatn, there have been no known catastrophic geological events that would necessarily have vastly reduced the long-run effective population size of small benthic and pelagic Arctic charr in Lake Thingvallavatn during the 3500–4000 generation scope of the sympatric phases of our historical models. Thus, it is also conceivable that population sizes have been very large, and such a scenario is very interesting, if indeed the contemporary population sizes do roughly reflect historical effective population sizes. In the pure sympatric scenario, population sizes much greater than 104.5 are inconsistent with the observed level of genetic differentiation (Figure 4). With such large, but nonetheless realistic, population sizes, an allopatric phase in the evolution of these morphs has to be invoked because otherwise drift is ineffective at generating the observed level of genetic differentiation.
The nature of the presumed barrier to gene flow between small benthic and pelagic Arctic charr is unknown. However, given the suggestion that this barrier might be very strong, it is tempting to speculate. Given the apparent opportunity for interbreeding, the maximum levels of gene flow between the populations that are consistent with the observed population-genetic statistics, that is, m<10−4 (Figure 4) are very low. Considering the near-universal tendency of male salmonids to sneak fertilizations (for example, Sigurjónsdóttir and Gunnarsson, 1989; Foote et al., 1997; Blanchfield et al., 2003) and further that this behavior is probably rather indiscriminate (Garcia-Vazquez et al., 2002), we should consider these rates of gene flow very small indeed. In situ observations of spawning behavior of the small benthic and pelagic Arctic charr morphs in Lake Thingvallavatn have not been made, although we do know that both male and female large benthic Arctic charr behave agonistically toward small benthic charr during spawning activities (Sigurjónsdóttir and Gunnarsson, 1989), and mate choice experiments have not been conducted. Hybrids of small benthic and pelagic charr from Lake Thingvallavatn have intermediate head morphologies (Skúlason et al., 1989a; Eiríksson, 1999), but much of the range of phenotypic consequences of hybridization remains to be investigated. We speculate that strong natural selection against intermediate morphologies and behaviors could be the barrier to gene flow.
Our historical models of the demography of the small benthic and pelagic Arctic charr morphs in Lake Thingvallavatn represent major simplifications of the potential complexity in the processes that generated contemporary patterns of genetic variation. We view our treatment of the problem as a generation of expected patterns of genetic variation along two informative slices through a vast parameter space. Indeed, variation along axes we have not investigated has probably been important. For example, we did not investigate variation in the size of the ancestral population nor did we investigate variation in the number of generations that have passed since the colonization (or colonizations) of Lake Thingvallavatn by Arctic charr. Our estimates of the numbers of generations is likely an overestimate, given the geological time frame (Sæmundsson, 1992) and the present life histories (Jonsson et al., 1988) of small benthic and pelagic Arctic charr in Lake Thingvallavatn. With respect to our ultimate interpretation that the observed genetic differentiation is consistent with high isolation, an overestimate of the number of generations that have passed renders this result, if anything, conservative. Given these unexplored axes of variation, we must qualify our finding that observed patterns in the distribution of genetic diversity among the morphs is consistent with very restricted gene flow, as just that, that is, consistent; we have not shown conclusively that gene flow has been low. Approaches do exist that may provide further insight into the joint effects of multiple variables (for example, the sizes of ancestral populations and the date of divergence), such as likelihood-based coalescent techniques (Hey and Nielsen, 2007) and approximate Bayesian computation (Beaumont et al., 2002). However, the information content of non-sequence data for this kind of inference is generally low, and it is unclear that the signatures of multiple historical processes in allele frequency data can necessarily be separated analytically.
At the larger spatial scale, the primary limitation of our analysis is the limited availability of the reference populations. Ideally, the analyses we undertook would have been conducted with multiple candidate ancestral populations for each population of the focal, that is, small benthic morph. We obtained such samples wherever possible. However, in many cases, no such (known) populations exist. More broadly, our definition of reference populations is based largely on necessity, in particular on the availablity of anadromous or otherwise unspecialized geographically proximate populations. We do not know from what morphs or specific populations any contemporary Arctic charr populations in Iceland were derived, and furthermore, such populations do not necessarily still exist. In particular, all morphs of Arctic charr that currently exist in Lake Thingvallavatn seem highly specialised and hence probably quite derived, especially relative to the probable ultimate ancestral form (likely an anadromous fish capable of colonizing Iceland from a glacial refugium that was probably located in Europe; Brunner et al., 2001). These limitations arising from incomplete knowledge and availability of candidate ancestral populations certainly hinder interpretations of the origins of any particular population of any particular morph. However, this level of uncertainty probably does not greatly affect our general conclusion that genetic variation is largely partitioned geographically, but not according to current and recent hydrological arrangements, and consequently that much of the evolution of the small benthic morph must have occurred multiple times.
Connectivity among Arctic charr populations in Iceland is generally low. Consequently, the evolution, and indeed the repeated evolution, of the small benthic morph has probably been relatively unconstrained by a homogenizing effect of gene flow. As a consequence of very low or non-existent contemporary and recent patterns of gene flow, genetic variation in Icelandic Arctic charr populations is structured much more by Cartesian overland distance than by contemporary waterway distance. This suggests that the genetic effects of demographic process during the colonization of Iceland still prevail in the patterns of genetic variation within Icelandic Arctic charr. Combined with the lack of association between morph and patterns of genetic differentiation, this suggests that repeated adaptive evolution of the small benthic morph is much more plausible than the alternate extreme hypothesis that the form evolved once, and is now widespread because of subsequent dispersal. Even within Lake Thingvallavatn, the level of genetic differentiation between the most genetically similar morphs, that is, the small benthic and pelagic forms, is consistent with low levels of gene flow, as a consequence of the very low rate of drift that will occur in very large populations. Relative to the range of degrees to which gene flow can act as an evolutionary constraint (for example, Jain and Bradshaw, 1966; Stearns and Sage, 1980; Hendry et al., 2002) and given the extensive ecological opportunity for diversification in a recently deglaciated landscape (Schluter, 2000), it seems clear that both selective and non-selective evolutionary forces have been important determinants of contemporary patterns of diversity in Icelandic Arctic charr.
We thank the farmer at Mjoanes, Johann Jonsson, for all his help with sampling from Thingvallavatn as well as the numerous field assistants who helped with the sampling of small benthic and reference Arctic charr populations. Xia Yue and Sheena Morrissey provided technical assistance with laboratory methods. Funding was covered by a RANN'IS grant to SSS, a RANN'IS grant from the Icelandic Research Fund for Graduate Students and a William and Ann Brock doctoral scholarship to BKK, an NSERC Discovery Grant to MMF and by NSERC PDF and PGS D awards to MBM.
The authors declare no conflict of interest.