Ecotype Differentiation in the Face of Gene Flow within the Diving Beetle Agabus bipustulatus (Linnaeus, 1767) in Northern Scandinavia

The repeated occurrence of habitat-specific polyphyletic evolved ecotypes throughout the ranges of widely distributed species implies that multiple, independent and parallel selection events have taken place. Ecological transitions across altitudinal gradients over short geographical distances are often associated with variation in habitat-related fitness, these patterns suggest the action of strong selective forces. Genetic markers will therefore contribute differently to differences between ecotypes in local hybrid zones. Here we have studied the adaptive divergence between ecotypes of the water beetle Agabus bipustulatus along several parallel altitudinal gradients in northern Scandinavia. This water beetle is well known for its remarkable morphological variation associated with mountain regions throughout the western Palaearctic. Two morphological ecotypes are recognised: a montane type with reduced flight muscles and a lowland type with fully developed muscles. Using a multilocus survey of allozyme variation and a morphological analysis with landmark-based morphometrics, across thirty-three populations and seven altitudinal gradients, we studied the local adaptive process of gene flow and selection in detail. Populations were sampled at three different elevations: below, at and above the tree line. The results indicate that the levels of divergence observed between ecotypes in morphology and allele frequencies at α-Glycerophosphate dehydrogenase relative to those shown by neutral molecular markers reflects local diversifying selection in situ. Four main lines of evidence are shown here: (1) A repeated morphological pattern of differentiation is observed across all altitudinal transects, with high reclassification probabilities. (2) Allele and genotype frequencies at the α-Gpdh locus are strongly correlated with altitude, in sharp contrast to the presumable neutral markers. (3) Genetic differentiation is two to three times higher among populations across the tree line than among populations at or below. (4) Genetic differentiation between ecotypes within independent mountain areas is reflected by different sets of allozymes.


Introduction
Gene flow is traditionally seen as the key factor that holds gene pools of local populations together, homogenising the genetic variation of interbreeding populations and opposing the effects of drift and local selection [1][2]. Recent empirical evidence challenges the universality of this view; site-specific selection and non-random mating have been shown to cause strong differentiation and/or reproductive isolation over small geographical scales [3][4][5][6][7]. As a result, focus is shifting towards understanding the evolutionary mechanisms driving adaptive versus neutral genetic differentiation within species where ecotypes exhibit sharp phenotypic differences across environmental boundaries over distances well within the dispersal capacities of a species [8][9][10][11][12][13]. The repeated occurrence of habitat-specific ecotypes within a widespread species is now taken as evidence for adaptation and, where gene flow can occur between ecotypes, that the process of divergence is driven by strong natural selection [5,[14][15][16]. Ultimately this knowledge will help to resolve one of the more contentious issues in evolutionary biology on how species or ecotypes may evolve and diverge in the absence of geographical barriers, namely, sympatric speciation. To date few credible examples of this process have been demonstrated although the number is growing [11,[16][17][18].
In general, ideal systems for exploring the potential role of diversifying selection for ecotypic divergence are found where common environmental factors, varying across some gradient, repeatedly lead to the existence of similar phenotypic forms in independent geographic localities [3,[12][13]. In this study, we focus on the water beetle Agabus bipustulatus (Linnaeus, 1767) one of the most common water beetles in the western Palearctic. It shows a remarkable morphological and geographical variation strongly associated with high altitude mountain regions [19][20][21][22][23]. Two morphologically distinct ecotypes are recognised [19,24] (Fig. 1): a montane type with reduced flight muscles and a lowland type with fully developed flight muscles [25][26][27][28]. The lowland form is known to be a strong flier able to migrate over the distances typically separating populations below the tree line, above the tree line and between mountains see [26][27][29][30]. In contrast, the reduced flight muscles of the montane form severely limit its dispersal capacity [28]. This form may migrate passively via downstream drift between populations within the same watercourse but migration between mountain tops is highly unlikely in the absence of flight. In addition to these disjunct distributions of morphological forms/ecotypes, the a-Glycerophosphate dehydrogenase locus (a-Gpdh), a major regulating enzyme in the metabolic pathway generating energy for flight [31][32][33], has been shown to display strong divergence between populations in the lowland areas and the alpine zones above the tree line. The locus has thus been identified as a candidate locus subject to selection mediated by factor(s) related to altitude and climate [21]. A phylogeographic analysis of the A. bipustulatus complex also confirmed that the associations between the a-Gpdh genotypes and altitude are observed across the entire west Palaearctic region [22] and that the montane ecotype has polyphyletic origins across the region [23].
These repeated and parallel distributions of morphological and enzyme phenotypes make A. bipustulatus an ideal system to study local environmental adaptation in the absence of gene barriers [12]. The Scandes mountain range provides an ecological theatre where these processes can be studied in more detail and on a scale where populations are more connected. Intermixed populations containing a few specimens of either ecotype have occasionally been observed above and below the tree line [28]. Looking across the Scandes mountain range, the distribution pattern of these ecotypes is, however, striking -populations on mountain peaks consist primarily of the montane form whereas all populations below the peaks consist almost entirely of the lowland type [19].
The purpose of this study was to analyse whether habitatspecific selection in A. bipustulatus is sufficiently strong to impede local gene flow across several independent parallel environmental gradients in Northern Scandinavia. Specifically, we analysed patterns of variation in a-Gpdh and morphological traits as well as in putatively neutral allozyme markers to assess genetic differentiation among populations across one longitudinal and seven altitudinal environmental gradients. These gradients span over three distinct climatic zones ranging from boreal through subalpine to alpine.
If in situ divergent selection is important for creating and maintaining ecotypic differentiation in the presence of countervailing gene flow in A. bipustulatus, we expect (i) sharp changes in morphological and a-Gpdh phenotypes, i.e. the traits subject to selection, across all tree lines. Since selectively neutral loci should show random patterns of variation among populations, especially given the potential for migration among populations, then if selection is imposing reproductive isolation (ii) the strongest differentiation in neutral markers should be observed between populations above and below the tree line both within and between mountains. Furthermore, individuals of different ecotypes within mountains should be less differentiated than individuals of the same ecotype from different mountains. (iii) If this is occurring independently in each location, the combination of neutral loci reflecting this differentiation should be specific to each mountain site.

The water beetle
Agabus bipustulatus is a dark medium-sized water beetle with a Palearctic distribution. The lowland ecotype is found in many types of habitats. In the boreal region, it prefers smaller streams, and springs and pools with little vegetation and stony bottoms. Above the tree line, the montane ecotype is most common in shallow lakes with no fish, chiefly with stony bottoms and some marginal vegetation [20].
The timing of male accessory gland enlargement indicates that the mating season is between May and September [34]. Eggs are laid either within submerged vegetation or in soft dry ground. Females generally carry one spermatophore, although two have been observed [35]. Females can lay between 600-800 eggs from June to August and have been observed to lay fertile eggs 12 months after the last copulation [34][35]. The eggs mature from early August to late April after which they hatch into first instar larvae. Agabus bipustulatus has three larval instars [36][37], all of which can be found almost year round in the lowlands of northern Sweden. In late autumn it is rare to find any first or second instar larvae above the tree line, where the third instar prevails [28]. After the third and final instar, the larvae enter the pupal stage which lasts for 42 days (at 11.3uC) [34]. Both adults and larvae are able to overwinter [20].
Water beetles have two major means of dispersal. For shortdistance dispersal (i.e. within a pond/lake), they swim or, if dispersal occurs along running water, passively drift with the current. The second mean of dispersal, flying, is mainly used for longer distances (i.e. between ponds/lakes). Main period for longdistance dispersal occurs in late summer [20]. The capacity for flight is not equally important for all species of water beetles and they show variation in flight muscle development, from absent to fully developed [25]. The mechanism behind this within-species variation has not yet been resolved. One plausible reason could be that a partial to full histolysis of the flight muscle occurs when flight capacity is obsolete. However, a study by Jackson [25] shows that no beetles with abnormal flight muscles could be induced to fly and the development of flight muscles and their support appeared to be arrested at an early stage. Her study also indicated a positive relationship between flightless forms and colder environments. Although as far as we know no study has documented the genetic control of flight muscle development. Studies of dispersal distances in A. bipustulatus are scarce but have shown that lowland forms with fully developed flight muscles disperse over distances up to several hundred meters [30]. Low mark-recapture numbers in areas with ponds separated by less than a couple of hundred meters indicate that the beetles can perform longer flights [34,38]. Other studies of similar sized water beetles have shown that they readily disperse several kilometres and distances of up to 25 km have been recorded [29]. Bilton [39] gives evidence for long-distance flying while Jackson [25] has tied flying ability to the development of the flight muscles. Low recapture rates with the studies cited above strengthen the idea that flight in the lowland form of A. bipustulatus is migratory as suggested by Southwood [40] and Dingle [41].

Sample collection
A total of 1279 individuals from 33 populations were sampled within seven altitudinal transects and one reference latitudinal transect in the Swedish Lapland during September and early October, 1998 (Fig. 2, Table 1). This sampling time is just after peak migration, which maximises the probability of finding both types of beetles in all habitats. Adult individuals were collected with a water bag net near rocks or overhanging vegetation at several locations along the margins of large water bodies, or throughout the entire water area in small water bodies. The specimens collected for the enzyme assays were kept alive and transported in damp moss to the laboratory where they were then separated by sex and frozen at 270uC. Specimens used in the morphological analyses were kept in 95% ethanol in a refrigerator at 4uC. Populations used in the morphological analyses are marked in italics in table 1.

Sampling design
The single latitudinal transect begins c. 100 km to the east of the mountain region and runs eastward towards the Baltic Sea coast for c. 200 km. This lowland region (LR) transect contains three sampled areas (named A-C), of which A has three populations and B and C have two each (Fig. 2). The seven altitudinal transects within the mountain region (MR) are located within different mountain districts, hereafter named D-J. All districts except D are separated by at least 15 km and have a mean diameter of 1-6 km. Area D has a diameter of 15 km but does not overlap with any other area. Districts I and J share the same mountain valley. Each of the altitudinal transects are divided into three elevation zones for sample collecting; below the tree line (hereafter referred to MB), at the tree line (MT) and above the tree line (MA). The tree line is defined as the highest continuous altitudinal level of Betula pubescens ssp. czerepanovii. It constitutes a major ecosystem shift between boreal and alpine zones [42][43]. All populations sampled above the tree line are from large shallow tarns, except Djuptjärn (I3), which is a deep water-turbine shaft. Population sites at the tree line vary from small cold spring outlets (Lillbraket-XIII, J1) and larger tarns (i.e. Kråkberget-IV-V (E1-4) and Unna Suojal, D3) to a large drilled shaft hole (Rundtjä rn, G1). The habitats below the tree line vary extensively, but most represent stable water bodies. The small number of specimens and the apparent instability of small shallow bodies of water with a high risk of drying-out, i.e. the habitats at Kråkberget-IV (E1), Vä stansjö-XVII (G3), Gröndal-XXI (H2) and Atoklinten-III (H6), suggest that these populations were newly founded.

Morphological analyses
Geometric morphometrics utilise the spatial covariation of homologous landmarks to calculate affine (partial warps) and nonaffine (uniform) shape components. These components can be used to calculate differences between morphology within and between ecotypes or the associations between shape and other parameters (e.g. altitude or genetic variation at specific loci (here a-Gpdh)) [44][45][46]. The shape component is a decomposed description of one specimen in relation to a consensus shape of all specimens in the analysis. Each landmark (geometric position) of a specimen is scaled (centroid size = 1), rotated, translated and aligned to minimise differences between samples (230 specimens in total, Table 1). The Procrustes metric (a) was set to null [47], which gives equal weights to partial warps at all spatial scales. The uniform component estimated as described by Bookstein [48] was included in the weight matrix.
We used eight landmarks, representing the beetles' right side ( These landmarks were selected to provide coverage of the shape variation observed in pronotal length, width and its anteriolateral variation relative to the elytra between the ecotypes [19,49]. Landmarks were not used on the head as its position relative to the rest of the body may obscure the analysis [50][51]. A Summa Sketch III (Summagraphics) graphics tablet was used to capture the geometric position of each landmark relative to a Cartesian coordinate. To assess the repeatability and minimise measurement errors of the geometric morphometric analyses, each specimen was measured three independent times as described by Lessells and Boag [52].
Morphological variation within the total sample was analysed with a Principle Component Analysis (PCA) of the partial warps and the uniform components in the Relative warps program v 1.20 [53]. The obtained components (relative warps) are only used to describe the shape variation and were not included in any statistical analysis [45][46]. The relationship between morphology (partial warp score and CS) and the number of copies of the a-Gpdh 100 allele and altitude in each sex were analysed using a Partial least square (PLS) analyses. This multivariate modelling analysis describes the relationship between sets of dependent (partial warp scores) and predictor variables (genetic variation, altitude and CS). The analyses were carried out separately for the two sexes in the SIMCA-S version 6.01 program [54].
Distribution pattern of the ecotypes were calculated by the posterior probability, based on the overall shape (partial warp scores), that a selected specimen of an ecotype belongs to a particular group (sample localities/altitudinal level). We thus use the over-all shape variation to estimate the number of montane ecotypes present at the different altitude levels rather than a visual determination. This analysis is done with a discriminate function analysis (DFA) and subsequent classification analysis. We expect to find specimens that are ''misplaced'' between ecotypes along the altitudinal gradient as also observed by Eriksson [28].

Genetic analyses
Prior to analyses of population structure, we tested whether the assayed loci could be assumed to behave in a selectively neutral manner. Neutrality can be assumed if the confidence intervals of F IP (the correlation of genes within Individuals relative to Population) overlap or lie between the confidence intervals of F IP calculated over loci [57]. The over-population confidence intervals of F IP were estimated from 5000 bootstrap randomisations. Per locus confidence intervals of F IP were calculated from jackknife re-sampling. Genetic diversity of each locus over population was calculated using Nei's unbiased estimator [58] in the FSTAT v2.93 software [59]. In order to access the genetic differentiation of neutral markers between populations above and below the tree line, both within and between mountains, we utilise a three-level hierarchical analyses of gene frequencies. The strength of the population differentiation was estimated via pairwise F ST calculations between populations at different elevation levels in mountain regions.
The relative levels of differentiation within and among altitudinal gradients in the mountain region (MR, areas D-J) were estimated from three-level hierarchical analyses of gene frequencies [60]. The following variance components were calculated from a three-level nested ANOVA: total variance s 2 t À Á , among areas s 2 a À Á , among populations within areas (s 2 b ), among individuals within populations s 2 c À Á , and within individuals (error term, s 2 w ). The estimators of genetic differentiation were obtained from these components as: F IP is the correlation of genes within Individuals relative to Population, F PA the correlation of genes within Population relative to Areas, F AT the correlation of genes within Areas relative to Total, F IA the correlation of genes within Individuals relative to Areas, F PT the correlation of genes within Populations relative to Total, and F IT the correlation of genes within Individuals relative to Total. The analysis was carried out using SAS v8.1 software with the ''NESTED'' procedure. Ratios of sums were used to obtain per locus and overall estimations [61]. Both the Mdh and Mnr loci were excluded from this analysis since scores from the Mdh locus in Braket-XIV and from the Mnr locus in Gieravardo and Kråkberget-IV were not considered to be reliable. To examine levels of genetic differentiation between elevations in the mountain region, population samples were designated as above (MA), at (MT) and below (MB) the tree line (Table 1). Genetic differentiation among populations within and between areas were estimated with pairwise F ST at their respective elevations; calculated groups are MAxMA, MTxMT, MBxMB, MAxMT, MBxMT, and MAxMB. Population differentiation is expected to be higher among populations above the tree line (MAxMA) than among populations below and at the tree line (MTxMT, MBxMB), since the populations of flying forms should be more connected than the populations of non-flying forms. For similar reasons, genetic differentiation between MBxMT elevation levels should be lower than between MAxMT and MAxMB within mountains. Estimates were carried out according to Weir and Cockerham [61] and analysed in the FSTAT v2.93 software [59].
To estimate the impact of individual loci on the genetic differentiation between subsamples of ecotypes on either side of the tree line population structure was estimated with multivariate ordination of subpopulations by principle component analysis (PCA) of allele frequencies with the PCAGEN software [62]. This analysis is independent of the assumption of Hardy-Weinberg equilibrium. From the PCA analysis, a two-dimensional canonical plot of the first two principal components was produced. The PCA result is used to discriminate between site-specific (nonallopatric) and allopatric origin of divergence.This was done on the total material (not including the a-Gpdh loci) and independently for the mountain areas F to J where the population sample-size made the comparison statistically meaningful ( Table 1). The purpose of this latter analysis was to test the prediction that if the two ecotypes have evolved in situ at each site, then the actual alleles contributing to the differentiation between ecotypes should be specific to each tree-line.

Ecotype distribution
The shape variation described by the first five relative warps in both sexes visualised the shape difference between A. bipustulatus lowland and mountain ecotypes (Fig. 1) and combined explained 85.3% and 84.5% of the total variation for females and males, respectively. Significant morphological discrimination (partial warp scores) between ecotypes was detected within the mountain region elevation levels (Males: Wilk's lambda = 0.009, p,0.001; Females: Wilk's lambda = 0.008, p,0.001). The morphological variation was strongly influenced by elevation level, which explained 89% of the total variation in the females and 88% in the males. If all single specimens were reclassified according to elevation level and significant discriminate function 83% of the females and 71% of the males were reassigned to their correct elevation group. Both ecotypes showed a similar pattern of reassignment (Fig. 3). The lowland ecotype specimens were ''misplaced'' between elevation levels (MA or MB) approximately three to four times more often than the montane ecotype.

Variation at the a-Gpdh locus
All loci surveyed, except a-Gpdh locus, had a none -significant overall population confidence interval (C.I.) of F IP that overlapped the C. I. over loci. The a-Gpdh locus has two alleles named 100 and 108. The a-Gpdh 100/100 genotype is significantly more frequent in populations above than at or below the tree line in the altitudinal gradient (Kruskal-Wallis test; H = 8.78, P = 0.032). The frequency of the a-Gpdh 108/100 genotype, on the other hand, decreases successively from below, to at and above the tree line (Fig. 4). No significant differences were observed between all three genotypes below the tree line (MB) and populations in the lowland region (LR). The unbiased genetic diversity estimator (H) of the a-Gpdh locus can be roughly divided into three categories: 1) genetic diversity within populations below the tree line (MB) and in the lowland (LR) varies between H = 0.4-0.5; 2) populations at the tree line (MT) H = 0.2-0.3 and 3) populations above the tree line (MA) H,0.1. However, exceptions exist; Rörmyrberget-III in the lowland has a genetic diversity of H = 0.32, and a few populations at the tree line in area G, I and J have a genetic diversity of 0.39 to 0.50, and no differences were observed across the tree line in area F and from the tree line to above in area E. Genotypic frequencies for each sampled population is given table 1.

Covariation of morphology and a-Gpdh across the tree line
The relationship between morphological variation (partial warp scores), altitude, CS (size) and the number of copies of the a-Gpdh 100 allele per population, were estimated with the Partial least square (PLS) analysis for both sexes independently. The analyses resulted in one significant (R1) component for each sex, with a correlation of 0.63 for females and 0.52 for males (Fig. 5). The significant components explained 17.6% of the morphological variation and 39.8% of the genetic variation in the a-Gpdh locus for females and 14.8% and 26.9% respectively for males. The predictive power (explained variation) of the PLS models for the Figure 3. Location of Agabus bipustulatus ecotypes in relation to sampled area and elevation in northern Sweden. Morphological data were generated from spatial covariation of homologous landmarks that calculate affine (partial warps) and nonaffine (uniform) shape components as described by Bookstein [48]. Dispersal patterns implied via a reclassification model following a discriminate function analysis (DFA) of the shape parameters are as follows:

Genetic structure
The genetic diversity over all populations was normally distributed and varied from 0.153 to 0.404 with an overall mean and standard deviation (SD) of 0.2616.07. Below the tree line, genetic diversity was similar in both the LR and MR. At the tree line, genetic diversity was similar in both the LR and MB regions. The diversity increased slightly at the tree line and then decreased moderately above it but these differences are not significant (Kruskal-Wallis test; H = 5.384, P = 0.146).
Significant positive deviations in genetic expectations under random mating were observed within 55% of the populations, as seen in the F IP values per population ( Table 1). The F IP values were normally distributed and varied between 20.032-0.449 over all populations (between 20.032-0.407 in the lowland region and 0.059-0.449 in the mountain region), with a mean and standard deviation of 0.20360.12.

Local genetic differentiation
In the mountain region, genetic differentiation between populations within areas and between populations within the region showed a moderate level in the hierarchical analysis, F PA = 0.080 and F PT = 0.075, respectively (Table 2). Most genetic variance is found among individuals within populations, s 2 c~3 8:9%, and 10.1% was found among populations within areas s 2 b À Á : Genetic differentiation was not detected between areas (F AT = 0.0 and s 2 a~0 , respectively). Pairwise F ST calculations, not including the a-Gpdh locus, were used to estimate genetic differentiation within and between different mountain areas, and between the three elevation classes: below (MB), at (MT) and above (MA) the tree line. These calculations showed that the mean pairwise F ST values between the elevation classes are higher between than within mountain areas. For elevation classes within mountain areas, the genetic differentiation increases from F ST = 0.026, between MB and MT, to F ST = 0.071 between MT and MA ( Table 3). The increase is less dramatic when elevation classes were compared between mountain areas, MBxMT, F ST = 0.051 and MTxMA, F ST = 0.078. A similar pattern of differentiation is also seen when elevation classes MB and MA are compared within and between areas, F ST = 0.065 and 0.075, respectively. Genetic differentiation among populations below and at the tree line is similar within mountains in the mountain region, whereas the genetic differentiation is always higher among populations when the montane ecotype above the tree line is included in comparisons, thus F ST = 0.050 within and 0.062 between mountains.
These results are congruent with the Principal components analysis (PCA) of the population structure, which explained 53% of the total genetic variation (not including the a-Gpdh locus) with a global F ST value of 0.097 (Fig. 6). Allele frequencies along the first significant component separated the two ecotypes with an overlap of the populations at the tree line, explained 36% of the genetic variation. The second axis was significant and explained 17% of the variation. The gradual genetic change is very strikingly but some population stands out, e.g. Gä utavardo (F3) and Kråkberget-VI (E3), as being more similar to the genetic composition of their ecotype counterpart.
The five loci with highest impacts (loading scores) separating the two ecotypes in the PCA analyses of mountains F to J are, in descending order: Area F: Ak, Mdh, Mnr, IDH, Est-2; G: IDH, Est-1, Est-2, Me, Mdh; H: Ak, IDH, Mdh, Mnr, Est-2; I: Ak, HK, Me, Mdh, Mnr, and area J: Ak, IDH, Est-1, Mdh, Mnr. Clearly, the alleles and loci that reflected ecotype separation in any one mountain area differed from the alleles and loci that separated the ecotypes in the other areas. This is strong support for in situ evolution of ecotypes within each mountain area.

Discussion
Here we will argue that the divergence between A. bipustulatus lowland and montane ecotypes in morphology and a-Gpdh genotype reflects independent local adaptive differentiation in situ across altitudinal gradients via four main lines of evidence: (1) Repeated patterns of morphological differentiation are observed in all independent gradients across the tree line, with high reclassification probabilities, at elevation levels (MA, MT and MB). Minor admixing of ecotypes occurs at all altitudinal levels. (2) The allelic and genotypic variation at the a-Gpdh locus is strongly correlated with altitude and habitats across the tree line, a pattern not observed in molecular markers deemed as neutral. (3) Molecular markers indicate that the strongest genetic differentiation is observed between ecotypes and that this differentiation is more pronounced between mountain areas than within, implying that gene flow may be impeded as a result of local adaptation. (4) The genetic differentiation between ecotypes, estimated from the enzyme systems, is reflected through different sets of loci in each independent mountain area. Evidently, the observed ecological transitions within A. bipustulatus across altitudinal gradients take place over short distances and are likely to be associated with variation in habitat-related fitness. Drotz et al. [23] have documented these associations between flightlessness and a-Gpdh variation in A. bipustulatus across the entire western Palaearctic region. The smaller scale process documented here could thus shed new light on the high morphological and taxonomic diversity within this species that has been discussed for more than 120 years [19,21,28].

Flight muscle hydrolysis and dispersal
The interpretation of the morphological and genetic adaptation across the tree line will depend on whether the dispersal potential of the ecotypes exceeds the realised gene flow [10,13]. Moreover the dispersal potential of the ecotypes will influence the cline width and the strength of the selective forces needed to maintain differentiation across the local hybrid zone [63].
Similar to the results of Eriksson [28], who showed that only 5.6% of 356 specimens of the mountain ecotype were found below the tree line in Finnish Lapland whereas 94% of the specimens above the tree line had reduced flight muscles, the morphological reclassification model, in this study, found higher numbers of ''misplaced'' specimens between below and at the tree line than above it (Fig. 3). This indicates that movement of beetles within valley floors are more common than between mountain tops of different Valley floors, while the downward dispersal of the montane ecotype is about 3-4 times lower than the upward migration of the lowland ecotype.
A comparison of the upward and downward migration may shed some light on when the flight muscle histolysis occurs. If muscle autolysis, which is a common phenomenon within  Coleoptera [64], occurs prior to oviposition or in the pupal stage, dispersal is more or less geographically restricted within water systems [63]. This may be true for A. bipustulatus. If hydrolysis begins after oviposition we would observe more montane type specimens below the tree lines in the mountain region, since more specimens should be able to fly. Our results indicate that the hydrolysis starts prior or during the pupal stage. Dispersal of the mountain ecotype should be more or less restricted to short dispersal events within lakes and/or downstream movements within water systems. These assumptions are supported by the morphological data for e.g. Atoklinten-II (H5), Atoklinten-IV (H7) and Braket-XIV (J2) which consist of individuals of extreme mountain type, high number of copies of the a-Gpdh 100 allele and come from different mountain ridges above the tree line (Fig. 5). These observations are interesting, since Atoklinten-I-III (H3, H5-6) is genetically closer to each other than to Atoklinten-V (H4) and Atoklinten-IV (H7) which belong to different water systems. This genetic similarity is also seen between Braket-XIV (J2) and Lillbraket (J1) which are only 1 km apart on the same mountain slope and share the same water system (Fig. 6). In contrast, to Djuptä rn (I3) and Kråkberget-VII (E4) which originate from the same altitude, but consist of less extreme specimens of the mountain type or are intermixed with the lowland type. They have a lower frequency of the a-Gpdh 100 allele. Gä utavardo (F3) is more similar to the lowland type morphological and genetically than to the mountain type (Fig. 6). Dispersal from the valley floor upwards across the tree line is therefore assumed to be mainly done by active flight.

Genetic structure
Large significantly positive F IP values were observed in most samples of both ecotypes of A. bipustulatus. These values are within the range observed in other water beetles, e.g. [22,56,65]. A reasonable explanation for this is common short-range dispersal behaviour of water beetles. Using marked beetles, both Süselbeck [34] and Davy-Bowker [38] showed that the autumn increases and spring decreases in population sizes are at least in part due to adults emerging from the pupae stage, immigration and emigration, respectively. They also noted that the beetles often only made one or few long distance migration bouts during the early stage of their life cycle after which they tend to stay in one specific region within a lake or smaller pond [29,34]. Specimens of A. bipustulatus tend to aggregate near large stones, where shore vegetation creates cover, or under large rocks near the shore line. However, beetles can be found all over the water body and bottom areas in smaller lakes. There is therefore a higher probability to sample different ''family'' aggregates due to admixture in larger habitats. This effect might explain the higher mean population F IP values [66] above the tree line in relation to other altitudinal levels and the lowland region, since the habitats above the tree line are larger than at or below it. Only two out of 13 smaller populations with fewer than 15 specimens had a significant F IP deviation (Table 1).

Adaptation across the tree line
The distance between populations below, at and above the tree line within each mountain area is short and well within the dispersal range of water beetles [29]. In spite of this, large genetic and morphological differences were detected between populations. Different neutral enzyme systems reflected these changes between ecotypes on the different mountains (Fig. 6). These observations indicate a constraint to gene flow and an independent adaptation of populations above the tree line (MA). This is supported by the variance components analysis ( Table 2) that showed that the genetic variance was higher within populations and mountain areas across the tree line than between collecting areas.
Levels of differentiation between populations in the lowland region (LR) and below the tree line (MB) in the mountain region (MR) are of similar magnitude implying similar levels of connectedness (Fig. 6). Gene flow within the mountain region, implied by pairwise F ST , takes place more easily between stable and temporary habitats along the Valley floor (MB) up to the tree line (MT) than between areas (Table 3). Genetic differentiation between populations across the tree line increases by a factor of three within areas and by 1.5 between areas in the mountain region. The estimated contribution of gene flow, not including the a-Gpdh locus, is therefore more pronouned from the Vally floor up past the tree line (Fig. 6). However, at the tree line populations tend also morphologically to be a mixture of specimens that are either of one or the other ecotype intermixed with few specimens of the other type (Fig. 5). The adaptive divergence consequently is most profound in populations above the tree line (MA). Gene flow does not, however, manage to conteract the effect of selection [8][9] and the morphological transition across the tree line must be abrupt, since few lowland ecotype specimens are found above the tree line (Fig. 3). The selection acting directly or indirectly on the a-Gpdh 108 allele must be strong [12,63,[67][68].
Studies of wing dimorphic species or comparisons between flying and non-flying species have often shown extensive genetic differentiation between populations or species. Within the flightless stream dwelling waterstrider Aquarius remigis, F ST values of 0.46 are reported among streams [69]. Lower mean heterozygosity values are reported for non-flying than flying waterstriders, H = 0.058 and 0.234, respectively [70]. Large F ST values have also been reported between flying and non-flying carabid species, F ST = 0.003-0.16 and 0.26-0.27, respectively [71][72]. A conclusion drawn from these studies is that population differentiation increases with reduced dispersal potential. Such a pattern supports our argument that flight muscle hydrolysis occurs prior to or during the pupal stage.
In conclusion, there is ample evidence that adaptive divergence maintains local ecotypic differentiation in spite of ongoing gene flow in A. bipustulatus. Thus the A. bipustulatus complex emerges as a promising example of non-allopatric evolution of in situ reproductive isolation [4][5][6][14][15][16]. In a future perspective we may ask the following question: is the morphological transition reversible between the parental generation and offspring? Is it possible to regulate/initiate the morphological transformation between ecotypes with temperature? Could this adaptive process described here explain the evolution of many high altitudinal water beetle species within the Palaearctic region [68]?