Morphological and genetic divergence in Swedish postglacial stickleback (Pungitius pungitius) populations
BMC Evolutionary Biology volume 11, Article number: 287 (2011)
An important objective of evolutionary biology is to understand the processes that govern phenotypic variation in natural populations. We assessed patterns of morphological and genetic divergence among coastal and inland lake populations of nine-spined stickleback in northern Sweden. Coastal populations are either from the Baltic coast (n = 5) or from nearby coastal lakes (n = 3) that became isolated from the Baltic Sea (< 100 years before present, ybp). Inland populations are from freshwater lakes that became isolated from the Baltic approximately 10,000 ybp; either single species lakes without predators (n = 5), or lakes with a recent history of predation (n = 5) from stocking of salmonid predators (~50 ybp).
Coastal populations showed little variation in 11 morphological traits and had longer spines per unit of body length than inland populations. Inland populations were larger, on average, and showed greater morphological variation than coastal populations. A principal component analysis (PCA) across all populations revealed two major morphological axes related to spine length (PC1, 47.7% variation) and body size (PC2, 32.9% variation). Analysis of PCA scores showed marked similarity in coastal (Baltic coast and coastal lake) populations. PCA scores indicate that inland populations with predators have higher within-group variance in spine length and lower within-group variance in body size than inland populations without predators. Estimates of within-group PST (a proxy for QST) from PCA scores are similar to estimates of FST for coastal lake populations but PST > FST for Baltic coast populations. PST > FST for PC1 and PC2 for inland predator and inland no predator populations, with the exception that PST < FST for body size in inland populations lacking predators.
Baltic coast and coastal lake populations show little morphological and genetic variation within and between groups suggesting that these populations experience similar ecological conditions and that time since isolation of coastal lakes has been insufficient to demonstrate divergent morphology in coastal lake populations. Inland populations, on the other hand, showed much greater morphological and genetic variation characteristic of long periods of isolation. Inland populations from lakes without predators generally have larger body size, and smaller spine length relative to body size, suggesting systematic reduction in spine length. In contrast, inland populations with predators exhibit a wider range of spine lengths relative to body size suggesting that this trait is responding to local predation pressure differently among these populations. Taken together the results suggest that predation plays a role in shaping morphological variation among isolated inland populations. However, we cannot rule out that a causal relationship between predation versus other genetic and environmental influences on phenotypic variation not measured in this study exists, and this warrants further investigation.
Morphological variation is a ubiquitous phenomenon and it exists at many taxonomic levels such as between species, within species and among individuals in populations. Depending on its relationship to the background environment, morphological variation may be either adaptive or non-adaptive [1, 2]. The degree to which the environment, phenotypic plasticity, and the underlying genetic architecture interact may also fundamentally influence morphological variation [3–5]. Moreover, divergent morphology is often associated with speciation events suggesting that phenotypic variation plays a role in the speciation process [6, 7]. Thus, a fundamental goal of evolutionary biology is to understand what governs patterns of morphological variation in natural populations.
In this study we explore patterns of morphological divergence and the genetic signature of recent isolation in populations of nine-spined stickleback (Pungitius pungitius) in northern Sweden. The nine-spined stickleback is a small euryhaline fish that inhabits a variety of freshwater and brackish environments in the temperate and subarctic regions of the northern hemisphere [8, 9]. In Sweden, nine-spined sticklebacks are widely distributed in the Baltic sea (3.5 - 8.8 psu ) and inhabit a variety of freshwater ponds, lakes and rivers .
Throughout its range, the nine-spined stickleback displays remarkable morphological variation, particularly in isolated freshwater populations where morphology appears to be shaped by a combination of isolation, drift and natural selection [11–13]. For example, sticklebacks from landlocked populations are generally larger than their coastal conspecifics . Also, loss of the pelvic girdle and dorsal, pelvic, and anal spines has been observed in several isolated populations [8, 9, 11–16], making it an excellent model for studies of parallel evolution.
The armour and spines of sticklebacks (Gasterosteidae) are used primarily as a defence against piscivorous fishes and many populations retain both characters in the presence of these predators . However there may be a cost in maintaining such defences in populations lacking predatory fishes, as reduced armour and spines may aid escape from invertebrate predators [18, 19]. If so, reduction or loss of these traits should be favoured in populations where invertebrate predators predominate [16, 19, 20]. Reversal of the loss of armour has also been shown in sticklebacks, suggesting that these defensive traits may be rapidly regained in populations based on the strength of natural selection experienced .
The topology of Scandinavia has been strongly influenced by recent climactic events related to its most recent deglaciation after the last glacial maximum approximately 11-10,000 calibrated years before present (ybp) [22, 23]. Prior to this time, the region was completely covered by an ice sheet. Fish populations originating from the Baltic Sea presumably became trapped in lakes and ponds formed by land uplift soon after the ice-melt [24, 25]. In northern Sweden, nine-spined sticklebacks occur in both single species lakes and lakes with mixed fish communities . Recently, predatory fish such as brown trout (Salmo trutta) and rainbow trout (Onchorynchus mykiss) have been introduced in a few lakes for recreational angling purposes, offering an ideal platform to explore the role of predation on morphological variation in nine-spined sticklebacks.
There are several goals to the current study. First, we investigate whether populations of nine-spined sticklebacks show genetic signatures of isolation in inland lakes and ponds isolated since the last glacial retreat. Here we investigate patterns of genetic differentiation, genetic variation, and bottlenecking using nine variable microsatellite loci. The second goal of the study is to investigate morphological divergence between populations that share common ecological conditions. To this end, we group populations according to whether they originate from the Baltic coast, coastal lakes, or inland lakes and ponds with or without predators. Finally, we explore the potential for natural selection to affect morphological evolution by comparing phenotypic divergence with neutral genetic divergence using a PST - FST approach; an alternative to QST - FST based on phenotypic variation in natural populations.
Summary statistics of microsatellite data per population are listed in Table 2. Fisher's exact tests indicated no significant departures from Hardy-Weinberg equilibrium for any locus in any population after application of a Bonferroni adjustment to correct for multiple tests . Tests of genotypic disequilibrium between loci were also non-significant after Bonferroni adjustment, supporting independent assortment of microsatellite loci. Analysis using MICROCHECKER detected possible null alleles at STN19 in GB and RT populations and potential null alleles or incorrect scoring due to stutter in STN196 in the RT population. Analysis of loci over all populations with LOSITAN suggests that loci are putatively neutral except locus STN196, which may be under positive selection, and STN19, which may be under balancing selection. We reran the LOSITAN analysis excluding inland populations and these results suggest that all loci in ancestral coastal populations are neutral.
A neighbour-joining phylogenetic tree using DA distances based on all microsatellite data demonstrates that Baltic coast and coastal lakes cluster together and are genetically more similar to one another than inland populations with or without predators (Figure 2a). Inland populations showed greater genetic distances, consistent with longer periods of isolation from coastal populations and each other (Figure 2a). Bootstrap support for the branching order of the inland populations is poor, which is in line with our expectation that these populations all originated within a short period of time after the last glacial retreat.
Population pairwise FST estimates demonstrate that coastal populations show lower levels of genetic differentiation than inland populations. Among coastal populations, pairwise FST values were low (range = -0.005 - 0.026, Table 3) indicating weak population structuring among these populations. In contrast, inland populations had high levels of genetic differentiation as evidenced by high pairwise FST values (range = 0.081-0.677). Pairwise standardized FST values gave similar results to pairwise FST values (Table 3). A nested within-group comparison between Baltic coast and coastal lake populations showed no significant differences in FST estimates (P = 1.00). A nested within-group comparison also shows no significant difference between FST estimates between inland lakes with or without the presence of predators (P = 0.57). However a significant difference between pooled coastal (coastal lakes, Baltic coast) and pooled inland (predator, no-predator) FST estimates is apparent confirming the visual impression from Figure 2 that inland populations show higher levels of within-population genetic differentiation (P = 0.001).
Inland populations had lower levels of within-population genetic variation than coastal populations. In coastal populations, all nine microsatellite loci were polymorphic, with 3-8 alleles segregating per locus. Baltic coast and coastal lake populations showed similar levels of observed heterozygosity (Ho), allelic richness (A), and gene diversity (hs) (P = 0.88, P = 0.89, P = 0.97 respectively, Table 4). Inland populations, on the other hand, had a maximum of 5 alleles segregating per locus and high levels of fixation ranging from 0-3 loci fixed for a certain allele per population. No two inland populations showed fixation for the same set of loci. Inland populations had significantly lower Ho, A and HS than their coastal counterparts (P = 0.01 for all variables, Table 4) but whether or not an inland population had predators present did not have a significant influence on either Ho, A, or HS (P = 0.31, P = 0.76, P = 0.47 respectively, Table 4).
We found evidence of recent population-level bottlenecking in two of the 18 populations investigated. Populations GB and VST both showed a significant heterozygosity deficiency based on our BOTTLENECK model criteria (Wilcoxon sign rank one tailed test: GB, P = 0.0137; VST, P = 0.0137). In all other populations there was no evidence for a departure from mutation-drift equilibrium implying that recent genetic bottlenecking is unlikely to account for reduced levels of genetic variation in inland populations.
Results of the analyses of the 11 morphological characters are summarized by population in Table 5 (See Methods for complete descriptions of morphological characters). Mean estimates of body length (SL) were highest for the inland populations without predators, but did not differ significantly from the other groups (Figure 3). Relative measures of head depth (per unit of standard length, SL-1) demonstrate that inland populations (with and without predators) have the most divergent morphology while coastal populations (Baltic coast, coastal lakes) show intermediate values for this trait (Figure 3). Relative head-eye length, on the other hand, show Baltic coast and inland lake without predator populations to have the most divergent morphology, with coastal lakes and inland lake - predator populations having intermediate values (Figure 3). Relative estimates of spine length (pelvic, anal, anterior dorsal, middle dorsal, posterior dorsal, girdle length and girdle width) consistently demonstrate significantly longer spines (SL-1) in coastal populations than in inland populations (Figure 3). The one measured meristic character, the number of dorsal spines, was variable in all populations with a mean of 9-10 spines (range 8-11) for all populations except BN, which had a mean of five spines (range 1-10) (Table 5). However, this character was not significantly different between groups (Figure 3).
A cluster analysis based on population means across all 11 morphological characters demonstrated strikingly similar morphology between Baltic coast and coastal lake populations (Figure 2b). Inland populations showed greater divergence in morphology than coastal populations (Figure 2b). Although there is some morphological clustering among inland populations that share similar predation regimes, there is also clustering of populations that have disparate predation regimes, suggesting that overall phenotypic appearance may evolve independently in these populations (Figure 2b).
Results of the principal components analysis (PCA) showed that the first two principal components had Eigenvalues greater than one and explained a combined 80.6% of the total variation of the full data. Therefore only these two axes were used in subsequent analyses. Principal components axis 1 (PC1) represented 47.7% of the variance and had high loading scores (> 0.40) for pelvic spine length, anal spine length, and all dorsal spine lengths after varimax rotation (Table 6). Principal components axis 2 (PC2) represented 32.9% of the variance and had high loading scores for standard length, head depth, head length, and pelvic girdle width (Table 6). Therefore, PC1 corresponds to a "spine length" function while PC2 corresponds to a "body size" function. Only pelvic girdle length had high loading scores for both PC1 and PC2 (Table 6).
Visualization of the PCA results reveals tight clustering of Baltic coast and coastal lakes in PC space (Figure 4), confirming their morphological similarity. Inland populations, on the other hand, showed much greater variation than coastal populations. Inland populations without predators showed greater variation in PC2 scores than in PC1 scores and the opposite is true for inland predator populations (Figure 3). Thus, the introduction of predators in inland populations appears to have reduced variation in body size, accompanied by an increase in the variability of spine length (Figure 4).
To test whether variation in PC scores differs significantly between groups and populations nested within groups, we constructed a nested multivariate analysis of variance (MANOVA) on PC1 and PC2. Our results show strong support for morphological divergence among groups and significant differences between populations nested within groups (Table 7).
PST - FSTcomparisons within groups
Within-group estimates of neutral genetic divergence (FST) and phenotypic divergence (PST) in Baltic coast and coastal lake groups demonstrated low estimates for each measure (Table 4, Figure 5). The Baltic coast group had PST > FST with non-overlapping 95% confidence intervals for both PC1 and PC2 suggesting that divergent selection on both spine length and body size is occurring in populations of this group. Coastal lake PST estimates, on the other hand, had overlapping 95% confidence intervals with FST (Table 4, Figure 5) which suggests that within this group, one cannot distinguish between the effects of natural selection and genetic drift on phenotypic appearance. Within the inland no predator group, we found P ST < FST for PC1 indicating stabilizing selection on spine length and the opposite trend, namely PST > FST for PC2, suggesting divergent selection in body size in this group (Table 4, Figure 5). Finally, inland populations with predators show PST > FST for PC1 and PC2 suggesting divergent selection on spine length and body size, but estimates of PST are more variable than any other estimate and confidence intervals overlap FST values, making a judgement as to whether natural selection or genetic drift is operating equivocal (Table 4, Figure 5).
To test whether our measures of PST are a robust and fair measure of QST, we investigated the estimate of the critical c/h2 ratio  for both PC1 and PC2 in cases where PST > FST and the inverse relationship (i.e. h2/c ratio) when PST < FST. Low values of these ratios indicate more robust support for PST to approximate QST in their various scenarios . If we assume that body size is moderately heritable so that h2 = 0.4 (see Discussion), and because c cannot exceed unity, using the lower 95% confidence limit for PST in cases where PST > FST yields low values of c/h2 critical values of 0.200 and 0.287 for Baltic coast populations PC1 and PC2 respectively (Table 4). This implies that a 70-80% lower c than h2 would be necessary to attribute the observed phenotypic divergence to drift. Thus, we may interpret our results to be robust under this scenario to the null expectation. Similarly, a low c/h2 critical value of 0.125 was found in inland no predator group for PC2 also demonstrating a robust signature of divergent selection within this group. Under the one instance where PST < FST in inland no predator PCI, the h2/c ratio was low at 0.283 and we may conclude that Q ST < F ST , i.e. a signal of convergent evolution/stabilizing selection on spine length in inland populations without predators.
In this study, we used molecular and morphological data to quantify differentiation between populations from coastal areas of the Baltic Sea and inland populations that have been isolated since the last deglaciation. In general, we find that Baltic coast and coastal lake populations are similar genetically and morphologically. Inland populations, on the other hand, show greater genetic and morphological divergence. The recent introduction of predators into some inland populations also appears to have altered the evolution of body shape and spine length as inland populations with predators varied considerably in spine length but showed much greater conservatism in body size than inland populations without predators. Taken together, these results suggest that nine-spined stickleback populations in northern Sweden are strongly influenced by a combination of recent glacial activity, isolation/drift, and natural selection.
Our results support the hypothesis that inland populations of nine-spined sticklebacks were isolated from coastal populations following the recent deglaciation of the Baltic Sea. Our first line of support for this hypothesis comes from the physical location of these sites which were among the first areas to be exposed after the retreat of the Scandinavian Ice Sheet circa 11 - 10 k cal. ybp [22, 23]. A plausible scenario is that these populations originated from Baltic stock and were isolated as the land uplifted from glacial rebound forming lakes and ponds [24, 25]. Additional support for this hypothesis is venerated by the observation that few natural populations of nine-spined sticklebacks are found in northern Sweden above or below the highest coastline exposed after the last glacial maximum .
In agreement with the hypothesis for isolation of inland populations, we found lower levels of genetic differentiation, as measured by Nei's DA, FST and standardized FST, in coastal sites and very high levels of genetic differentiation including the fixation of different alleles in many inland populations. Accordingly, we found higher levels of genetic variation as measured by HO, π, and hs in coastal populations than in inland populations. These last two findings are consistent with a pattern of isolation and genetic differentiation of populations through genetic drift and/or founder effects. Based on our results, these decreases in HO do not appear to be due to a recent bottleneck in inland populations, with the exception of VST where significant results of simulated genetic bottlenecking were detected. We also found that populations from coastal lakes appear to be genetically indistinguishable from those emanating from the Baltic Sea, suggesting a separation so recent that genetic differences have yet to become apparent. These results strongly suggest that coastal populations are closely related whereas inland populations have existed in isolation for quite some time, and that genetic drift is responsible for the divergence in allelic frequencies and fixation of certain alleles.
Our support for a reduction in genetic variability among inland populations echoes findings in freshwater populations of the closely related three-spined stickleback (Gasterosteus aculeatus) that generally have reduced genetic variability [29–31] as a result of isolation from presumably ancestral anadromous populations . Similar models for post-glacial colonization by nine-spined sticklebacks have also been proposed and confirmed for North American populations [8, 33]. A recent study of nine-spined sticklebacks by Shikano and colleagues demonstrated similar patterns of isolation from anadromous ancestral populations of this species in Europe . In their study, coastal populations from the Baltic sea had high levels of genetic variation as measured by allelic richness and heterozygosity when compared to populations from freshwater systems near the Baltic while higher levels of genetic differentiation, as measured by FST and DA, were prevalent in freshwater but not in coastal populations . The authors concluded that these patterns of genetic variation and genetic differentiation are consistent with postglacial recolonization of freshwater habitats, and subsequent isolation reducing variation in these populations through genetic drift and founder effects. Since Shikano and colleagues' study encompassed sampling from a much larger geographic area than studied here, it appears that this pattern of isolation is common in the Baltic region and that recent glacial history has greatly affected the current distribution of these fishes.
Our results show many morphological similarities between Baltic coast and coastal lake populations of nine-spined sticklebacks, mirroring the pattern of low levels of genetic differentiation. Given that coastal lakes have been isolated from the Baltic due to land uplift over a relatively short period of time (< 100 yrs), it is not surprising that there has been little differentiation in morphology compared to the relatively greater morphological divergence in inland populations. Moreover, piscivorous fish such as pike and perch are not currently detected in two of the three coastal lakes  but were likely present in these habitats prior to and during lake formation. Therefore morphological divergence in the absence of predation may be expected in these populations in the future, but perhaps not over the short time since isolation from the Baltic.
Inland populations of nine-spined sticklebacks displayed higher morphological diversity as compared to coastal populations. We also found evidence of morphological differences in both body size and spine length with respect to the presence of fish predators. Additionally, we found highly divergent morphological variation in some inland populations. For example, the populations with the two most divergent morphologies, BN and RT, show the greatest differences in mean body size. The BN population also has a reduction in the number of dorsal spines and dorsal spine length, and all but one individual completely lacked pelvic spines. This pattern of pelvic spine loss and an overall reduction in spines has been demonstrated in several species of stickleback [16, 19, 34] and may be due to an ecological escape from predation pressure, an ion deficiency related to calcification and bone deposition and/or increased invertebrate predation pressure in these populations [16, 20, 35].
Nine-spined sticklebacks that exist in ponds where they are the only fish species can obtain much larger body sizes than their coastal counterparts, presumably because of the absence of fish predation combined with interspecific competition for resources and/or fecundity selection [11, 12]. Confirming this prediction, we observed a trend for individuals hailing from inland populations without predators to be larger, on average, than both Baltic populations and inland lakes with predators. A similar pattern of larger size encountered in populations that either lack predators or have non-gape limited predators has been shown in three-spined sticklebacks  and brook sticklebacks (Culaea inconstans) , strongly suggesting that predation limits body size in this and other species of sticklebacks.
Despite having a smaller range in body size than their counterparts that hail from predation-free lakes, nine-spined sticklebacks from predator lakes have much greater variation in spine length. Although spine lengths are significantly smaller in inland populations compared to coastal populations, there is a trend for inland populations with predators to show longer spines (per unit of body length), on average, than inland populations without predators. The implication, therefore, is that spine length has increased in these populations as a result of the recent introduction of predatory fish. However the great variation in spine length exhibited among these populations within this group suggest that these populations may be responding to predation pressure differently. This great variation spine length could be caused by differing predation regimes within these lakes due to differences in predator communities, population densities of predators, differences in the relative exposure time to predators, and/or the recent extinction of predators within at least one lake. Although our results strongly suggest a role of predation, we cannot rule out that other ecological factors not measured in this study such as water chemistry may affect morphological variability in predator lakes.
PST - FSTcomparisons
Quantitative comparisons between morphological and neutral genetic divergence as estimated by QST - FST, have been used as metric to investigate the potential for natural selection to influence morphological variation in many populations or species [38, 39]. The analogous comparison based solely on phenotypic data (i.e., PST - FST) has been highly criticized, particularly because it is difficult to tease apart the variance in phenotype attributable to environmental or genetic effects in wild populations. Thus, some advocate that QST - FST comparisons should only be performed under controlled conditions in common garden experiments . While such approaches are preferable, we argue that there can be some value to PST - FST comparisons in natural populations. For example, in a recent meta-analysis that compared estimates of QST - FST from different types of studies, estimates of PST from wild populations do not yield higher estimates than studies that use either broad or narrow sense estimates of additive genetic variation . Secondly, studies of PST in wild populations show meaningful variance among populations where a common garden approach may not be easily applied [27, 29, 41]. For example, in order to quantify additive genetic variation in our study of 18 natural populations with all potential crosses and multiple family groups taken to the F2 generation would be a feat of herculean proportions in terms of time, scale and expense. Finally, in light of the criticism of PST as a substitute for QST, it should also be kept in mind that common garden estimates of QST may be inappropriate to compare to FST because the genetic basis of the phenotypes on which selection may potentially act may be partly genetic but non-additive (i.e. epigenetic), or environment-dependent. Thus we believe that PST - FST studies in natural populations do have some merit although we advocate caution in its interpretation.
Acknowledging the aforementioned concerns, we compared PST estimates for both spine length (PC1) and body size (PC2) within groups to estimate the relative influence of natural selection and genetic drift on the evolution of morphological phenotypic variation within our groups of interest. Despite the great morphological similarity within the Baltic coast group, it demonstrated PST > FST in both spine length and body size. Given the high gene flow experienced within this group, the most likely explanation is that these slight morphological differences could be explained by phenotypic plasticity in response to local environmental variation. However it should be noted that our FST estimates are so low that any variation in morphology would likely be greater than genetic divergence and so this result should be viewed cautiously. In coastal lakes, PST ≈ FST for these estimates making distinctions between selection and genetic drift equivocal. In inland lakes without predators, we found a robust FST - PST pattern of convergent evolution and stabilizing selection on reduced spine length strongly suggesting that reduced spine length is advantageous in single species lakes, potentially to aid escape of fish from invertebrate predators. A robust pattern of divergent selection and local adaptation on body size is also evident in inland no-predator lakes. Finally, great variation in body size and spine length within the inland predator group may also be indicative of divergent selection on these traits in these populations but confidence intervals are wide making a realistic determination difficult. Taken together, these results show opposite and robust patterns of convergent versus divergent selection on spine length based on the presence or absence of predation, strongly implicating the recent exposure to predators as a significant factor shaping phenotypic differences between these populations.
The degree to which our main focal phenotypic characters, body size and spine length, differ in their phenotypic response to selection is not currently known in nine-spined stickleback but may vary among different populations due to different environment conditions and standing levels of genetic variation [3, 5]. Studies on other fish species show that fish predators potentially can induce phenotypic changes in body shape and morphology without necessarily changing the background genetic structure [42–44]. For example, one study showed that cues of a predatory fish induced a deeper body and longer dorsal spines in a sunfish (Lepomis gibbosus) . The genetic component of the traits studied here remain unknown as previous studies investigating genetic differences among populations of nine-spined stickleback using a common-garden approach have not yet reported heritablities for these traits [13, 45]. Among three-spined sticklebacks, morphological traits of predator-naive and predator-sympatric populations demonstrate high values of heritability for body shape (h2 = 0.67, 0.94) and length (h2 = 0.92, 0.82) and moderate values for relative spine length (h2 = 0.34, 0.39) . Other studies of three-spined sticklebacks have shown that body size is moderately heritable (h2 = 0.42)  or that heritability may vary from negligible (h2 = 0.007) to moderate (h2 = 0.313) due to environmental factors such as rearing salinity . Taken together, these results suggest that in most cases morphological traits such as body size and spine length are moderately to highly heritable in three-spined sticklebacks, but that the ratio between genetic and environmental variance (and hence the heritability) may not be equal across populations.
Understanding the causes underlying morphological differences between populations is a fundamental question in evolutionary biology. Our power to explore differences in morphology among inland populations with disparate predation regimes is significantly enhanced by our knowledge of predator stocking and 'fisher knowledge' based on over 100 years of data collection within inland lakes. Our two findings that 1) fish hailing from inland no predator lakes appear to be convergent selection on small spine length and divergent selection on body size and 2) fish from predator lakes have a highly variable response in spine length and body size strongly suggests that the recent introduction of salmonids in these lakes have influenced morphology in nine-spined sticklebacks in these populations. However, our current study cannot determine whether additional factors such as water chemistry [20, 35] and invertebrate predation  can contribute to phenotypic variation in our populations of stickleback. Additionally, the amount of environmental variation and heritability of phenotypic traits within and between populations is outside of the scope of this study but certainly warrant additional investigation along these lines. This research also highlights the importance of studying morphological and genetic variation on a fine ecological scale in order to determine the environmental factors responsible for shaping the shared and unique features of evolutionary histories.
Field sampling of nine-spined stickleback populations occurred during May - August 2007 and June 2008 in the counties of Västerbotten and Västernorrland, northern Sweden (Table 1, Figure 1). In total, 18 locations were visited including eight populations from the Baltic Sea coast consisting of five populations from coastal bays (Baltic coast) and three populations from lakes situated within a km of the Baltic Sea (coastal lakes). An additional 10 populations were sampled from inland lakes and ponds that have been isolated from the Baltic Sea for more than 9000 cal. ybp (Figure 1). Of the 10 inland lakes and ponds, five are single species lakes where nine-spine sticklebacks are the only fish species recorded (inland - no predators, Table 1). Stickleback populations in the remaining five inland lakes and ponds have a known history of predation (inland - predators). Predators such as brown trout (Salmo trutta) and rainbow trout (Oncorhynchus mykiss) have been stocked in these lakes for recreational angling purposes. Records show that these lakes (except HM) have been stocked as early as 1959 and many have been stocked extensively in the 1980s and 1990s (Table 1). No records of the presence of these fish precede the known stocking dates. Records indicate that one of our sampling sites (SKT) was stocked with brown and rainbow trout in 1985, but surveys after 2000 indicate that these introduced species are no longer present in this lake. Nevertheless, this lake was classified as a predator lake based on the very recent exposure to predation.
The selection of locations was based on information in the PIKE Database . The database contains records of fish community composition and introductions for more than 19,000 Swedish lakes. Information about fish species is based on interviews with local fishermen and surveys using a minimum of eight minnow traps and/or four multi-mesh gill net sets (bottom set, monofilament nylon and measured 30 × 1.5 m with mesh size ranging between 5 and 55 mm) .
Fish collection and morphological measurements
Adult nine-spined sticklebacks were caught using baited minnow traps and landing nets. Fish were immediately placed on ice and transported to the laboratory for genetic and morphological analyses. There, specimens were defrosted and morphological traits were measured with a digital calliper. Traits measured included standard length (SL, tip of snout to tip of caudal peduncle), head depth (top of head to bottom of head through centre of eye perpendicular to SL), head-eye length (top of head to centre of eye perpendicular to SL), length of pelvic spine, length of anal spine, and length of anterior, middle and posterior dorsal spine (Figure 6a). The total number of dorsal spines was also counted. For individuals with even numbers of dorsal spines, middle dorsal spine length was an average of the length of the two middle spines. After measurements were taken, the right pelvic fin (or the caudal fin in the BN population) was clipped for genetic analyses, and individual fish were stained in order to measure the pelvic girdle.
To prepare the staining of the pelvic girdle, fish were first fixed according to the following scheme: 95% EtOH for 24 h, 70% EtOH for 24 h, 50% EtOH for 24 h, 20% EtOH for 24 h and 24 h in distilled H2O. Fish were then transferred to a 0.4 g/l aqueous solution of KOH containing 0.425 g of Alizarin red to stain bony parts. Fish were kept in the solution for 3 h and were then transferred to 0.4 g/l aqueous solution of KOH without dye. After 24 h, fish were placed in 50% iso-propanol for long-term storage. After staining, the length and width of the pelvic girdle was measured using callipers (Figure 6b).
Total genomic DNA was extracted from finclips of all individuals in each sampled location using a DNeasy ® Blood & Tissue Kit (product #69581, Quigen, Valencia, USA). Nine polymorphic microsatellite markers (Pbbe1125, Stn19, Stn49, Stn96, Stn148, Stn163, Stn173, Stn196, Stn198) previously used to characterize P. pungitius  were individually labelled with florescent primers. Microsatellite markers were amplified using polymerase chain reaction (PCR) in a 10 μl reaction using a multiplex PCR kit (Qiagen Inc.). PCR conditions for each reaction using the multiplex kit were 5 μl of Qiagen Master Mix (containing HotStarTaq DNA polymerase, dNTPs and 3 mM MgCl2), 0.625 μM of each primer, and 3 μl of genomic DNA. Temperature profile for thermal cycling: 96°C for 15 min followed by 35 cycles of 94°C for 30 sec, 60°C for 90 sec, 72°C for 90 sec, and a final 20 min extension at 72°C. Fragment lengths of PCR products were analyzed using CEQ™ 800 Genetic Analysis System (Beckman Coulter Inc., Brea, USA).
Microsatellite data were analyzed with ARLEQUIN v. 3.11  to test for Hardy-Weinberg equilibrium (Fisher's exact test) and for genotypic disequilibrium for pairs of loci within populations (Fisher's exact test). Loci were checked the scoring errors and the presence of null alleles using MICROCHECKER v 2.2.3 . We also tested all loci for neutrality using the program LOSITAN using a multi-step mutation model and 10,000 replicates of data collection [52, 53].
A single microsatellite locus, Stn49 showed a mutation consisting of a one base pair insertion. This mutation is common in eastern European lineages of P. pungitius  and it was included in all subsequent analyses except for analyses conducted with MICROCHECKER, LOSITAN, and BOTTLENECK because this locus does not follow a typical dinucleotide stepwise mutation pattern .
Nei's genetic distance (DA)  was used to construct a phylogeny based on microsatellite data for the 18 nine-spined stickleback populations using the program POPULATIONS v. 1.2.30 . Levels of bootstrap support were obtained from 10,000 replicates and the resulting neighbour-joining tree was drawn with POPULATIONS.
Genetic differentiation, as measured by FST , was calculated between each population pair with ARLEQUIN. Significance of global pairwise FST values was estimated with ARLEQUIN using 10,000 permutations. In order to compare measures of genetic differentiation corrected for differences between populations with different levels of allelic richness and heterozygosity, we calculated standardized FST values with GENODIVE v. 2.0b18 .
We tested for differences in within-group genetic variation (Baltic coast, coastal lakes, inland lakes with or without predators) and between coastal (pooled Baltic coast, coastal lakes) and inland lakes (pooled predators, no predators). Our estimates of genetic variation included observed heterozygosity (Ho), expected heterozygosity (He) allelic richness (A), and global FST values with FSTAT v 18.104.22.168  using the "comparison among groups of samples" subfunction. Significance was ascertained for these parameters over 10,000 permutations.
To investigate whether populations bear the signature of a recent fluctuation in population size, we investigated heterozygosity deficiency and heterozygosity excess with respect to gene diversity using the program BOTTLENECK v. 1.2.02 with 1,000 replications of a two-phase mutation model . The two-phase model of mutation consisted of mostly of one-step changes with a low percentage of multistep changes (9:1), recommended by Di Rienzo et al.  for microsatellites.
We first tested for differences between groups (Baltic coast, coastal lakes, inland no predators, inland predators) using a general linear mixed model (GLMM) using group as a fixed factor and population nested within group as a random factor to ascertain whether groups significantly differed in standard length and the number of dorsal spines. All remaining morphological measurements were first divided by standard length to obtain a relative morphological measurement per unit of standard length (SL-1) prior to GLMM analysis. A Tukey-Kramer honestly significant difference post-hoc analysis was conducted to determine which groups gave rise to significant differences.
A morphological cluster analysis was conducted on population means of all 11 morphological characters including the number of dorsal spines. Principal components analysis (PCA) was performed on 10 morphological characters (excluding the meristic character) and a varimax rotation  was applied to all components that had Eigenvalues greater than 1. A nested multivariate analysis of variance (MANOVA) was then used to test for differences in morphological PC axes between all groups, with populations nested within groups. Since many individuals in the BN population lacked dorsal, pelvic and anal spines or had highly reduced spine lengths and therefore may have contributed a disproportionate amount of variation to the PCA, we constructed the PCA and all subsequent analyses with and without the BN population. However, all significant differences in MANOVAs between groups were detected irrespective of whether BN was included, so we present only the results with BN included.
Calculation of PST
In can be shown that (for diploid species, assuming purely additive gene action) the quantity
(where V A, b and V A, w are the morphological additive genetic variance components between and within populations) is expected to take the same value as FST calculated from the genes that determine the trait [62, 63]. Hence, if FST is calculated from neutral markers such as the microsatellites we analyzed, then for a quantitative trait under divergent selection QST > FST and under stabilizing or convergent selection QST < FST. Unfortunately, the common garden experiments needed to estimate the additive genetic variance components of a character are difficult to realize for many species. Therefore, QST is often approximated by PST, which is calculated directly from the total phenotypic (i.e. combined genetic and environmental) variance components :
Obviously, how well PST approximates QST depends on whether the total phenotypic variance within and between populations reflects the additive genetic variance. The heritability h2 relates the additive variance within a population to the total variance: V A, w = h2 V w . We can similarly envision a "population level heritability" c so that V A, b = c V b . If we are willing to make the simplifying assumptions that both h2 and c are equal for all populations, then [27, 63]:
Algebraic rearrangement then yields the relation between QST and PST expressed as a simple function of the ratio c/h2:
Thus, we can calculate PST and evaluate which value the ratio c/h2should take for QST to equal FST, that is:
In other words, if for some trait we find that PST < FST, we can calculate which ratio c/h2 would generate this result when in fact QST = FST. In the reverse situation when PST > FST, we can take the inverse of the c/h2 ratio (i.e., h2/c) to test the robustness of our estimates. If the required ratio is extreme, i.e. c < < h2 we may be confident that the conclusion based on a PST - FST comparison would not be altered if we had estimated QST directly. Therefore, we report "uncorrected" values of PST in the results, and use the above formula when interpreting the robustness of our results.
To calculate P ST , we partitioned phenotypic variance in within- and between population components by ANOVA (assuming unequal variances) in Matlab (version 7.5, MathWorks, Natick, Massachusetts, USA). We then calculated the within-group mean pairwise PST of all pairs of populations. If there are n populations in the group, there are n (n-1)/2 pairs of populations. For example, for the n = 5 coastal bay populations, PST is calculated as the mean of 10 pairs of populations. For between-group PST we calculated the mean pairwise PST of all pairs of populations where each pair consists of a population from the one group and a population from the other group. If there are n1 populations in the one group and n2 populations in the other, there are n1 n2 between-group pairs. Because each population is part of several pairwise comparisons in the calculation of PST, we used bootstrapping to calculate 95% confidence intervals.
Unless otherwise noted, all statistical tests were performed with JMP IN™ statistical software package version 7.0 (SAS Institute Inc. Cary, North Carolina, USA). Means are reported throughout text as ± one standard error of the mean (SE).
Rundell RJ, Price TD: Adaptive radiation, nonadaptive radiation, ecological speciation and nonecological speciation. Trends in Ecology & Evolution. 2009, 24: 394-399. 10.1016/j.tree.2009.02.007.
Schluter D: The Ecology of Adaptive Radiation. 2000, Oxford, U.K.: Oxford University Press
Pfennig DW, Wund MA, Snell-Rood EC, Cruickshank T, Schlichting CD, Moczek AP: Phenotypic plasticity's impacts on diversification and speciation. Trends in Ecology & Evolution. 2010, 25: 459-467. 10.1016/j.tree.2010.05.006.
Hodgins-Davis A, Townsend JP: Evolving gene expression: from G to E to G × E. Trends in Ecology & Evolution. 2009, 24: 649-658. 10.1016/j.tree.2009.06.011.
McGuigan K, Nishimura N, Currey M, Hurwit D, Cresko WA: Cryptic genetic variation and body size evolution in threespine stickleback. Evolution. 2010, 65: 1203-1211.
Langerhans RB, Giford ME, Joseph EO: Ecological speciation in Gambusia fishes. Evolution. 2007, 61: 2056-2074. 10.1111/j.1558-5646.2007.00171.x.
Nosil P, Crespi BJ: Experimental evidence that predation promotes divergence in adaptive radiation. Proceedings of the National Academy of Sciences of the USA. 2006, 103: 9090-9095. 10.1073/pnas.0601575103.
Aldenhoven JT, Miller MA, Corneli PS, Shapiro MD: Phylogeography of ninespine sticklebacks (Pungitius pungitius) in North America: glacial refugia and the origins of adaptive traits. Molecular ecology. 2010, 19: 4061-4076. 10.1111/j.1365-294X.2010.04801.x.
Shikano T, Shimada Y, Herczeg G, Merilä J: History vs. habitat type: explaining the genetic structure of European nine-spined stickleback (Pungitius pungitius) populations. Molecular Ecology. 2010, 19: 1147-1161. 10.1111/j.1365-294X.2010.04553.x.
Samuelsson M: Interannual salinity variations in the Baltic Sea during the period 1954-1990. Continental Shelf Research. 1996, 16: 1463-1477. 10.1016/0278-4343(95)00082-8.
Herczeg G, Gonda A, Merilä J: Evolution of gigantism in nine-spined sticklebacks. Evolution. 2009, 63: 3190-3200. 10.1111/j.1558-5646.2009.00781.x.
Herczeg G, Gonda A, Merilä J: Rensch's rule inverted - female-driven gigantism in nine-spined stickleback Pungitius pungitius. Journal of Animal Ecology. 2010, 79: 581-588. 10.1111/j.1365-2656.2010.01665.x.
Herczeg G, Turtiainen M, Merilä J: Morphological divergence in North-European nine-spined sticklebacks (Pungitius pungitius): signatures of parallel evolution. Biological Journal of the Linnean Society. 2010, 101: 403-416. 10.1111/j.1095-8312.2010.01518.x.
Ziuganov VV, Zotin AA: Pelvic girdle polymorphism and reproductive barriers in the ninespine stickleback Pungitius pungitius (L.) from northwest Russia. Behaviour. 1995, 132: 1095-1105. 10.1163/156853995X00478.
Shapiro MD, Summers BR, Balbhadra S, Aldenhoven JT, Miller AL, Cunningham CB, Bell MA, Kingsley DM: The genetic architecture of skeletal convergence and sex determination in ninespine sticklebacks. Current Biology. 2009, 19: 1140-1145. 10.1016/j.cub.2009.05.029.
Blouw DM, Boyd GJ: Inheritance of reduction, loss, and asymmetry of the pelvis in Pungitius pungitius (ninespine stickleback). Heredity. 1991, 68: 33-42.
Wooton RJ: A Functional Biology of Sticklebacks. 1984, London: Croom Helm
Marchinko KB: Predation's role in repeated phenotypic and genetic divergence of armor in threespine stickleback. Evolution. 2008, 63: 127-138.
Reimchen TE: Spine deficiency and polymorphism in a population of Gasterosteus aculeatus: an adaptation to predators?. Canadian Journal Of Zoology. 1980, 58: 1232-1244. 10.1139/z80-173.
Bell MA, Orti G, Walker JA, Koenings JP: Evolution of pelvic reduction in threespine stickleback fish: a test of competing hypotheses. Evolution. 1993, 47: 906-914. 10.2307/2410193.
Kitano J, Bolnick DI, Beauchamp DA, Mazur MM, Mori S, Nakano T, Peichel CL: Reverse evolution of armor plates in the threespine stickleback. Current Biology. 2008, 18: 769-774. 10.1016/j.cub.2008.04.027.
Björck S: A review of the history of the Baltic Sea, 13.0-8.0 ka PB. Quaternary International. 1995, 27: 19-40.
Björck S: The late Quaternary development of the Baltic Sea basin. Assessment of climate change for the Baltic Sea Basin. Edited by: Team BA. 2008, Berlin Heidelberg: Springer-Verlag, 398-407.
Ekman S: Djurvärldens Utbredningshistoria. 1922, Stockholm: Bonniers Boktryckeri
Englund G, Johansson F, Olofsson P, Salonsaari J, Ohman J: Predation leads to assembly rules in fragmented fish communities. Ecology Letters. 2009, 12: 663-671. 10.1111/j.1461-0248.2009.01322.x.
Rice WR: Analyzing tables of statistical tests. Evolution. 1989, 43: 223-225. 10.2307/2409177.
Brommer JE: Wither PST? The approximation of QST by PST in evolutionary and conservation biology. Journal of Evolutionary Biology. 2011, 24: 1160-1168. 10.1111/j.1420-9101.2011.02268.x.
Raeymaekers JAM, Van Houdt JKJ, Larmuseau MHD, Geldof S, Volckaert FAM: Divergent seletion as revealed by PST and QTL-based FST in three-spined stickleback (Gasterosteus aculeatus) populations along a coastal-inland gradient. Molecular Ecology. 2007, 16: 891-905.
Withler RE, McPhail JD: Genetic variability in freshwater and anadromous sticklebacks (Gasterosteus aculeatus) of southern British Columbia. Canadian Journal Of Zoology. 1985, 63: 528-533. 10.1139/z85-078.
Rafinski J, Banbura J, Przybylski M: Genetic differentiation of freshwater and marine sticklebacks (Gasterosteus aculeatus) of eastern Europe. Zeitschrift Fur Zoologische Systematik Und Evolutionsforschung. 1989, 27: 33-43.
Bell MA, Foster SA: The evolutionary biology of the threespine stickleback. 1994, New York: Oxford University Press
McPhail JD: Geographic variation in North American ninespine sticklebacks, Pungitius pungitius. Journal of the Fisheries Research Board of Canada. 1963, 20: 27-44. 10.1139/f63-004.
Reist JD: Variation in frequencies of pelvic pehnotypes of the brook stickleback, Culaea inconstans, in Red water drainage, Alberta. Canadian Field-Naturalist. 1981, 95: 178-182.
Giles N: The possible role of environmental calcium levels during the evolution of phenotypic diversity in Outer Hebridean populations of the three spined stickleback, Gasterosteus aculeatus. Journal of Zoology. 1983, 199: 535-544.
Reimchen TE: Trout foraging failures and the evolution of body size in stickleback. Copeia. 1991, 36: 1098-1104.
Zimmerman MS: A field study of brook stickleback morphology: multiple predators and multiple traits. Canadian Journal of Zoology. 2007, 85: 250-260. 10.1139/Z07-003.
Merilä J, Crnokrak P: Comparison of genetic differentiation at marker loci and quantitative traits. Journal of Evolutionary Biology. 2001, 14: 892-903. 10.1046/j.1420-9101.2001.00348.x.
Leinonen T, O'Hara RB, Arias JMC, Merilä J: Comparative studies of quantitative trait and neutral marker divergenece: a meta-analysis. Journal of Evolutionary Biology. 2008, 21: 1-17.
Pujol B, Wilson AJ, Ross RIC, Pannell JR: Are QST - FST comparisons for natural populations meaningful?. Molecular Ecology. 2008, 17: 4782-4785. 10.1111/j.1365-294X.2008.03958.x.
Sæther SA, Fiske P, Kålås JA, Kuresoo A, Luigujõe L, Piertney SB, Sahlman T, Höglund J: Inferring local adaption from QST - FST comparisons: neutral genetic and quantitative trait variation in European populations of great snipe. Journal of Evolutionary Biology. 2007, 20: 1563-1576. 10.1111/j.1420-9101.2007.01328.x.
Vøllestad LA, Varreng K, Poleo ABS: Body depth variation in crucian carp Carassius carassius: an experimtnetal individual-based study. Ecology Of Freshwater Fish. 2004, 13: 197-202. 10.1111/j.1600-0633.2004.00048.x.
Januszkiewicz AJ, Robinson BW: Divergent walleye (Sander vitreus)-mediated inducible defenses in the centrarchid pumpkinseed sunfish (Lepomis gibbosus). Biological Journal of the Linnean Society. 2007, 90: 23-36.
Brönmark C, Miner JG: Predator-induced phenotypic change in body morphology in crucian carp. Science. 1992, 258: 1348-1350. 10.1126/science.258.5086.1348.
Herczeg G, Gonda A, Merilä J: Predation mediated population divergence in complex behaviour of nine-spined stickleback (Pungitius pungitius). Journal of Evolutionary Biology. 2009, 22: 544-552. 10.1111/j.1420-9101.2008.01674.x.
Dingemanse NJ, Van der Plas F, Wright J, Réale D, Scharama M, Roff DA, Van der Zee E, Barber I: Individual experience and evolutionary history of predation affect expression of heritable variation in fish personality and morphology. Proceedings Of The Royal Society B. 2009, 276: 1285-1293.
Alho JS, Leinonen T, Merilä J: Inheritance of vertebral number in the three-spined stickleback (Gasterosteus aculeatus). PLoS One. 2011, 6: e19579-10.1371/journal.pone.0019579.
Appelberg M: Swedish standard methods for sampling freshwater fish with multi-mesh gillnets. Fiskeriverket information (Swedish board of fisheries report). 2000, 2000: 1-
Mäkinen HS, Valimaki K, Merilä J: Cross-species amplification of microsatellite loci for nine-spined stickleback Pungitius pungitius. Annales Zoologici Fennici. 2007, 44: 218-224.
Schneider S, Roessli D, Excoffier L: ARLEQUIN v2.1: an integrated software for population genetics data analysis. 2000, Genetics and Biometry Laboratory, Department of Anthropology, University of Geneva
van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P: MICRO-CHECKER: software for identifying and correcting genotype errors in microsatellite data. Molecular Ecology Notes. 2004, 4: 535-538. 10.1111/j.1471-8286.2004.00684.x.
Antao T, Lopes A, Lopes RJ, Beja-Pereira A, Luikart G: LOSITAN: A work bench to detect molecular adaptation based on a FST outlier method. BMC Bioinformatics. 2008, 9: 323-10.1186/1471-2105-9-323.
Beaumont MA, Nichols RA: Evaluating loci for use in the genetic analysis of population structure. Proceedings Of The Royal Society B. 1996, 263: 1619-1626. 10.1098/rspb.1996.0237.
Nei M, Tajima F, Tateno Y: Accuracy of estimated phylogenetic trees from molecular data. Journal of Molecular Evolution. 1983, 19: 153-170. 10.1007/BF02300753.
Langella O: POPULATIONS, 1.2.30. 1999, [http://bioinformatics.org/~tryphon/populations/]
Weir BS, Cockerham CC: Estimating F-statistics for the analysis of population structure. Genetic Research. 1984, 74: 215-221.
Meirmans PG, Van Tienderen PH: GENOTYPE and GENODIVE: two programs for the analysis of genetic diversity of asexual organisms. Molecular Ecology Notes. 2004, 4: 792-794. 10.1111/j.1471-8286.2004.00770.x.
Goudet J: FSTAT v. 22.214.171.124. 2002, [http://www2.unil.ch/popgen/softwares/fstat.htm]
Piry S, Luikart JL, Cornuet JM: BOTTLENECK: a computer program for detecting recent reduction in the effective size using allele frequency data. Journal of Heredity. 1999, 90: 502-503. 10.1093/jhered/90.4.502.
Di Rienzo A, Peterson AC, Garza JC, Vlades AM, Slatkin M, Freimer NB: Mutational processes of simple-sequence repeat loci in human populations. Proceedings of the National Academy of Sciences of the USA. 1994, 91: 3166-3170. 10.1073/pnas.91.8.3166.
Tabachnick BG, Fidell LS: Using multivariate statistics. 2001, Boston, MA: Allyn & Bacon
Spitze K: Population structure in Daphnia obtusa: quantitative genetic and allozymic variation. Genetics. 1993, 135: 367-374.
Merilä J, Crnokrak P: Quantitative trait and allozyme divergence in the greenfinch (Carduelis chloris, Aves, Fringillidae). Biological Journal of the Linnean Society. 1997, 61: 243-266.
Leinonen T, Cano JM, Makinen H, Merilä J: Contrasting patterns of body shape and neutral genetic divergence in marine and lake populations of threespine sticklebacks. Journal of Evolutionary Biology. 2006, 19: 1803-1812. 10.1111/j.1420-9101.2006.01182.x.
We thank Fredrik Engdahl, Jenny Lund, Joakim Moberg, and Stina Gustavsson for assistance in the field. We also thank Nate Jue, Per Ingvarsson and Melanie Monroe for technical advice and assistance. The Swedish Geological Survey (SGS) kindly provided 9,000 and 10,000 calibrated ybp coastline maps. This study received financial support from the Department of Ecology and Environmental Science and Umeå University.
KBM conducted all molecular and morphological analyses and drafted the manuscript. DL collected samples and measured morphological traits. FJ and GE organized collections and supervised morphological analyses. FB conducted PST calculations and participated in the design and coordination of the study. All authors read and approved the final manuscript.
Authors’ original submitted files for images
About this article
Cite this article
Mobley, K.B., Lussetti, D., Johansson, F. et al. Morphological and genetic divergence in Swedish postglacial stickleback (Pungitius pungitius) populations. BMC Evol Biol 11, 287 (2011). https://doi.org/10.1186/1471-2148-11-287