Multiple Origins of Elytral Reticulation Modifications in the West Palearctic Agabus bipustulatus Complex (Coleoptera, Dytiscidae)

The Agabus bipustulatus complex includes one of Europe's most widely distributed and common diving beetles. This complex, which is known for its large morphological variation, has a complex demographic and altitudinal variation in elytral reticulation. The various depth of the reticulation imprint, both in smaller and larger meshes, results in both mat and shiny individuals, as well as intermediate forms. The West Palearctic lowland is inhabited by a sexually dimorphic form, with shiny males and mat females. In mountain regions, shiny individuals of both sexes are found intermixed with mat individuals or in pure populations in central and southern areas, whereas pure populations of mat individuals are exclusively found in the northern region at high altitude. Sexual selection is proposed as a driving force in shaping this variation. However, the occurrence of different types of reticulation in both sexes and disjunct geographical distribution patterns suggest an additional function of the reticulation. Here we investigate the phylogeographical history, genetic structure and reticulation variation of several named forms within the Agabus bipustulatus complex including A. nevadensis. The molecular analyses recognised several well-supported clades within the complex. Several of the named forms had two or more independent origins. Few south European populations were uniform in reticulation patterns, and the males were found to display large variation. Reticulation diversity and population genetic variability were clearly correlated to altitude, but no genetic differences were detected among populations with mixed or homogenous forms. Observed reduction in secondary reticulation in female and increased variance in male at high altitude in South Europe may be explained by the occurrence of an additional selective force, beside sexual selection. The combined effect of these selective processes is here demonstrated in an extreme case to generate isolation barriers between populations at high altitudes. Here we discuss this selective force in relation to thermal selection.


Introduction
The elytral surface of most diving beetles is covered by a reticulation pattern impressed on the elytra, often consisting of small and large meshes -the primary and secondary reticulation [1][2][3]. The primary reticulation consists of a more or less regular fine meshwork, whereas the secondary reticulation is coarser with larger meshes of a more irregular, often oblong, shape. Sexual dimorphism is commonly found, especially in geographically disjunct areas at high altitude or gradually between longitudinal areas in e.g. the West Palaearctic region. Traditionally, these reticulation patterns have taxonomically been recognised as a set of characters used to identify and separate species, subspecies and varieties within several species complexes [1,[4][5][6][7][8]. However, little effort has been taken to understand the evolutionary significance of the sexual and altitudinal variation in elytral microsculpture and reticulation patterns occurring within a single species complex. This view changed when Bergsten and co-workers [9] nicely demonstrated an ongoing arms-race including female elytral sculpture and the number and size of adhesive male pro-and mesotarsal setae in several dytiscine genera. Miller [10] and subsequent studies [11][12][13][14] strengthened these observations and argued with support from the general model of sexual conflict [15][16] that the first step of a female response to increased mating cost is the occurrence of an increased female reticulation compared to males with tarsal modifications. Accordingly, these findings provide an explanation of the sexual dimorphism and increased modification of female reticulation in dytiscids. However, the observation that both males and females in some water beetles display a large variation in primary and secondary reticulation indicates that the microsculpture most likely has evolved for some other reason. Moreover, the occurrence of similar reticulation patterns at high altitude in different mountain ranges across large geographical areas also indicates the presence of an additional selective agent or common evolutionary history [17][18][19]. Sexual selection could therefore be a selective mechanism that enhances or contributes secondary to the diversification of reticulation patterns at high altitude in some species [20][21][22].
Clearly, several underlying evolutionary factors need to be considered in order to understand why different aberrant reticulation patterns, such as a more strongly impressed male secondary reticulation and more shiny females, occur at high altitude in contrast to the more homogeneous reticulation pattern normally observed at lower altitudes [19,[23][24]. Here we will focus on the phylogeographic history of different reticulation patterns and genetic differences between populations consisting of various reticulation forms in relation to altitude observed in the Agabus bipustulatus (Linnaeus, 1767) species complex.
This species complex is known to display a pronounced individual shape variation in combination with an elytral reticulation that varies both between and within sexes, especially at higher altitudes [1][2][3][4][5][6][7][8]13,[25][26]. The pronounced morphological, geographical and altitudinal variation is reflected taxonomically in the 23 available junior synonyms of uncertain taxonomic status [8]. The montane form named Agabus solieri Aubé, 1837 has a body-shape, especially seen on the pronotum, which is more slender than in the widely distributed topotypic lowland form bipustulatus [25,27]. Most of the other named forms are montane, and found in populations imbedded in areas dominated by the solieri-like form. Main differences between the forms, which taxonomically often have been dealt with as subspecies or varities of either bipustulatus or solieri, are the shininess of the elytra and different patterns in their secondary reticulation (Fig. 1). In particular, the kiesenwetterii form stands out since both sexes have wider meshes in the secondary reticulation, and the elytra are more shining relative to the regular solieri form. The geographical distribution of these two forms is interesting since the kiesenwetterii form commonly is found at high altitude in central and south Europe imbedded in populations dominated by the duller solieri form, whereas solieri is the only form occurring in North Europe [1]. These forms are suggested to be subjected to sexual selection which is seen in populations where males evolve wider tarsi and attached suckers in the presence of modified reticulated females [13]. In other known forms, only the male displays reticulation variation, whereas the females are homogenous for a given pattern.
The Iberian endemic Agabus nevadensis Hå. Lindberg, 1939 is here of interest due to its geographically restricted high-altitude occurrence in the Sierra Nevada and morphological similarity to A. bipustulatus. This local endemic has a reticulation pattern identical to that of the kiesenwetterii form in both sexes [28]. Its species status has been questioned and it is commonly included within the A. bipustulatus-complex [7,26].

Sample Collection
A total of 717 specimens representing six of the named forms of A. bipustulatus were collected from 37 sample sites across its West Palearctic distribution area in order to capture as large genetic and morphological variation as possible (Table 1, Fig. 2). All specimens used in the population genetic analysis were also included in the morphological analysis of the secondary reticulation.
In France, 14 A. bipustulatus populations were sampled near Grenoble, Guillestre, and Montcenisio. The type locality of A. solieri lies within the Grenoble region, and Guignot [7] described A. solieri var. falcozi from Montcenisio. Seidlitz [29] described A. solieri var. kiesenwetterii from Illyria, Piemonte and the Pyrenees, and it is according to Guignot [7] now also found it at high altitudes in the French Alps. One Italian locality, from which A. bipustulatus dolomitanus Scholz, 1935 is known, was visited [26]. Five populations were sampled in the Spanish Pyreenes near Vall de Boí in Catalunya, including the water system from which A. solieri pyrenaeus was described [30]. A total of 413 A. nevadensis specimens from eight populations were sampled near its type locality in central Sierra Nevada [28] (Table 1, Fig. 2).
To analyse the evolutionary background of different reticulation forms we utilised two A. nevadensis specimens and 30 A. bipustulatus specimens representing the reticulation pattern of 16 bipustulatus, two kiesenwetterii, three falcozi, one dolomitanus, one pyrenaeus and seven solieri forms. Individual specimens were collected to represent different secondary reticulation forms from different geographical regions with the aim to maximise the haplotypic (mtDNA) variation. The optimal outgroup taxa for the A. tristis species group is A. bipustulatus and A. nevadensis, according to molecular analyses of the Agabini, species from the A. nebulosus group [31]. Two A. nebulosus (Forster, 1771) specimens were therefore included, together with one A. affinis (Paykull, 1798) and one A. guttatus (Paykull, 1798). The former represents a basal clade of the subgenus Gaurodytes, whereas the latter belongs to the sister group of the A. nebulosus plus A. tristis groups [31]. Seven specimens of the A. tristis group not belonging to the A. bipustulatus complex were also included in the phylogenetic analysis ( Collected specimens that were used in the population genetic and phylogenetic analyses were transported alive in damp moss to the laboratory, where they were sexed and frozen at 270uC. Morphologically analysed specimens were stored in 96% ethanol in a refrigerator at 4uC.

Reticulation Patterns
The reticulation patterns of the sampled forms are as follows; in North Europe, from the high mountains of Scandinavia to the lowland region below the Alps and the Pyrenees, the A. bipustulatus complex is represented by a sexually dimorphic form with a more or less constant elytral reticulation pattern. The males are shinier than the females and the secondary reticulation of both sexes consists of longitudinally stretched meshes. In the montane   kiesenwetterii form both sexes have wide meshes in the secondary reticulation and the elytra are shinier than in topotypic A. bipustulatus from Sweden. Guignot [7] observed that falcozi males have wider secondary reticulation and are not as shining as those of kiesenwetterii, whereas the females are similar to the normal bipustulatus form. Only the males of pyrenaeus are shiny and have a more isomorphic secondary reticulation, and females have been described as similar to those of bipustulatus [30]. The latter information is not fully correct, however, and we will here demonstrate that the pyrenaeus females also display a similar variation as the males. The dolomitanus form is similar to falcozi, but it has the pronotal shape of lowland bipustulatus [32]. Franciscolo [32] also concluded, without any quantitative analysis, that the variation in secondary reticulation in kiesenwetterii and dolomitanus is so similar that they have to be conspecific. Individual variation in the secondary elytral reticulation was classified after the total sample was screened for common patterns. Seven categories (A -G) ranging from tightly packed longitudinally elongated meshes to more or less isomorphic meshes were recognised (Fig. 1). Of these categories, the A, B, C and G patterns cover most of the total surface of the elytra ( Table 2). The remaining categories D, E and F describe intermixed patterns of longitudinally elongated and isomorphic meshes. The A reticulation pattern is found in the bipustulatus and solieri forms. The B pattern is found in kiesenwetterii, falcozi and within the dolomitanus form. Within pyrenaeus we recognised three reticulation patterns: C, F and G.
The Simpson index (D), which measures the probability that two individuals randomly selected from a sample will belong to the same reticulation category, was used to calculate reticulation diversity per population [33].

Genetic Analysis
A total of eleven enzyme systems coding for nine loci gave quantitatively reliable results: Hexokinase (Hk), Triosephosphate isomerase (Tpi), a-Glycerophosphate dehydrogenase (a-Gpdh), Isocitrate dehydrogenase (Idh), Malate dehydrogenase (Mdh), Malic enzyme (Me) and Esterase (Est-1 & Est-2). In addition, Alkaline phosphate (Alp) was screened in A. bipustulatus from France and Italy, and so were Phosphoglucoisomerase (Pgi), Phosphogluconate dehydrogenase (Pgd) and Xanthine dehydrogenase (Xdh) in A. bipustulatus and A. nevadensis from Spain. Staining recipes and gel buffers are modified from those of Shaw and Prasad [34]. The beetles were prepared and the alleles named by migration distance to the most common allele as described in Drotz et al. [27].
Population structure was calculated as described by Weir and Cockerham [35] and genotypic linkage disequilibrium was tested for between loci per population. To evaluate if reticulation patterns influence population differentiation as assumed by sexual conflict in combination with assortative mating [36][37], we tested for differences in gene diversity, observed heterozygosity and F ST among groups consisting of specimens within populations with mixed and homogeneous reticulation patterns. Populations were grouped according to their variability in secondary reticulation observed as described above. Significant levels were calculated via a permutation scheme of 5000 replicates where the total sample is allocated at random to groups. The P-value is based on a twotailed probability test of the proportion of the randomised data sets that gives larger mean values than the observed. All analyses were conducted in the FSTAT v2.91 software [38]. Pearson's correlation coefficient was estimated between altitude and mean heterozygosity and F IS .

DNA Amplification and Sequencing
DNA was extracted from one hind leg or the thoracic flight muscles of the frozen or alcohol preserved beetles, with a Qiagen DNeasy protocol for animal tissues. Partial cytochrome b (Cyt b) sequences were amplified with the two flanking primers CB-J-10933 and CB-N-11367, partial cytochrome c oxidase subunit I (COI) with CI-J-2195 and TL2-N-3014, and cytochrome oxidase II (COII) with C2-J-3279 and C2-N-3661 as described in Simon et al. [39] using the following polymerase chain reaction (PCR) program: denaturation 94uC (90 s), annealing at 50uC (30 s) and extension at 72uC (60 s). This cycle was repeated 30 times, followed by an extension period of 5 minutes. The PCR products obtained were single clear bands with no signs of non-specific amplification. The amplified product was approximately 430 bp in length for Cyt b, 730 bp for COI, and 400 bp for COII. The product was run on a 1.0 % agarose gel and then removed from the gel and purified with a Jetsorb DNA extraction kit. Sequencing reactions were performed with the DYEnamic ET terminator kit. Each sequence was sequenced from both its 39 and 59 ends. Corresponding accession number for the three partial mtDNA genes are AF439354-55, AF439359, AF439362, AF439365, AF439367-68, AF439372-77, and AY535285-400.

Phylogenetic Analysis
The mtDNA sequences were aligned with the ClustalW multiple alignment option in BioEdit, version 4.8.10 [40]. No gaps were inserted within the alignment. Congruence among data sets was tested with the partition homogeneity test [41]. Unweighted parsimony analysis was performed by applying heuristic search with tree bisection-reconnection branch swapping. A total of 3000 searches were done with 100 replicates and ten random-addition sequence iterations per search started from a random tree. Branches were collapsed if branch length was minimum zero. All characters were nonadditive, and uninformative characters were excluded before the analysis. To evaluate if the data set is subjected to 'long branch attraction' we compared the strict consensus tree topology between two phylogenetic analysis; the first including all sequences from all outgroups and the second analysis including only the sequences from the ingroup as described by Bergsten [42]. The second dataset was also analysed with a maximum likelihood analysis. The Akaike Information Criterion in Modeltest [43] was used to find, the best evolutionary model given the data. Cladograms were rooted between the ingroup and outgroup according to Nixon and Carpenter [44]. Nodal support within the phylogenetic trees within the parsimony analysis was assessed with bootstrap percentage after 1000 re-sampling steps [45] and within the maximum likelihood analysis with jackknife percentage after 1000 re-sampling steps with 30% character deletion [46]. The number of extra steps required to collapse each clade was also calculated as described by Bremer [47]. Measures of how well all individual character fit on a phylogenetic tree is measured by the consistency index (CI). The average value is calculated by dividing the minimum possible number of steps by the observed number of steps on the tree. The retention index (RI) measures the amount of Table 2. Agabus bipustulatus classification of secondary elytral reticulation after the total sample was screened for common patterns.

Code Description
A Longitudinally elongated meshes more or less straight, narrow and a majority are long. Isomorphic meshes can be observed at the elytral apex.
B Longitudinally elongated meshes more or less straight, wide and a majority are long. Isomorphic meshes are observed commonly at the elytral apex and occasionally along the anterior parts of the elytra suture.
C Longitudinally elongated meshes more or less straight, wide and a majority are short.
D Longitudinally elongated meshes more or less straight, wide and a majority are long. Isomorphic meshes are observed from elytral apex to K or L of its length.
E Longitudinally elongated meshes more or less straight, wide and a majority are long. Isomorphic meshes are observed from elytral apex to K or L of its length and at the anterior parts of the suture.
F Longitudinally elongated meshes more or less straight to K of elytral length after that the meshes bend out towards suture J from apex. Isomorphic meshes at the elytral apex and along anterior parts of the suture.
G Isomorphic meshes cover the whole elytra, giving an impression that the longitudinally elongated meshes bend out towards suture J from apex.
Seven categories (A-G) ranging from tightly packed longitudinally elongated meshes to more or less isomorphic meshes were recognised in Agabus bipustulatus. Type A is found in the bipustulatus and solieri forms, type B in kiesenwetterii, falcozi and dolomitanus, and type C, F and G in pyrenaeus. synapomorphy expected from a data set that is retained as synapomorphy on a tree. The above analyses were run in both PAUP v.4.0b10 [48] and WinClada ver. 0.9.99 [49]. Tajima's Neutrality test for the mtDNA sequence data [50] was performed with DNAsp [51] and the assumption of a local or global molecular clock was tested with PHYHYP v. 1 [52]. The Kishino and Hasegawa [53] test was used to test for two competing trees of equal length. To evaluate differences between tree topologies of unequal length we performed a permutation test where character states are randomly reassigned within characters of all taxa. Length differences between all permutated trees generate a null distribution that is used to compare the observed length differences of the original topologies. Significant level is gained in both analyses from a two-tailed probability test in PAUP v.4.0b10 [48].
Reticulation patterns were not optimised as a character on any of the parsimonious fundamental trees. Instead we here utilized the strict consensus to study the within and between clade distribution of reticulation patterns over all parsimonious fundamental trees. This method is well suited to distinguish if reticulation patterns have evolved multiple times within different monophyletic clades [23].

Reticulation Diversity
The secondary reticulation varies noticeably within almost all populations from France, Italy and Spain (Table 3). In three populations (BipFra4, 6 and 7), all at an altitude below 2000 m a.s.l near Grenoble, all specimens of the same sex had similar reticulation patterns. The other populations, sampled at higher altitudes, consisted of a mix of A and B patterns in various frequencies with one or several individuals representing the D, E or F patterns. The exceptional Spanish population BipEsp5 did not have any males of neither the A nor B type, and instead we here found the C, F and G patterns dominating in both sexes (Table 3). Some geographically very close populations, separated by altitude, as BipFra5 and BipFra6 or BipFra13 and BipFra14, had different dominating male elytral patterns (A and B, respectively). In males the variation is significantly correlated with altitude (p,0.001, adj-R 2 = 60.2%, b = 0.79). Increased diversity of reticulation patterns is seen in several independently parallel areas in Spain and France (Fig. 3). In females, a similar pattern, however not significant, of high diversity at high altitude is also observed.

Genetic Variability
The within population variation in reticulation pattern made it difficult to assign a population into a given form. Populations were therefore divided into four groups based on the most frequent reticulation pattern in order to test for genetic differentiation between groups: (1) type A in all females and males (BipFra4, 6 and 7); (2) type A in all females and .50% type A in males (BipFra3, 8, 9 and 14); (3) type A in all females and .50% type B in males (BipFra9, 10, 12 and 13); (4) both sexes with variable reticulation patterns (BipFra1, 2, 5 and 11).
Agabus nevadensis displays a significant deviating pattern in observed heterozygosity (p = 0.0054), genetic diversity (p = 0.0008) and allele frequencies in relation to A. bipustulatus. The mean heterozygosity of A. nevadensis is H = 0.13460.044 and the mean number of alleles per locus is 1.9. Of the compared enzyme systems, the allele frequencies of Hk and Pgd stand out. The Pgd

Sequence Variation
Combined, the three mtDNA data sets COI, COII and Cyt b included 1236 base pairs. A total of 214 sites were parsimonyinformative over all taxa, 150 within the A. tristis group, and 33 within the A. bipustulatus complex. Of the 214 informative characters, 174 were from the third, 34 from the first, and 6 from the second position. Translation of nucleic acid to amino acid sequences revealed very few parsimony-informative characters; i.e. most of the variation seen represents silent mutations. The incongruence length difference test was not significant between any combination of the three data sets. Most of the character conflicts in the data are within the separate sets, and only 2.33-4.19 % increased homoplasy was found in the total data set including all taxa.
Tajima's test of the neutral mutation hypothesis of the A. nevadensis and the A. bipustulatus complex sequences showed no significant deviation (p.0.10) between the amount of segregating/ polymorphic sites and the average number of nucleotide differences, number of segregating sites (S) = 52, Theta per sequence = 13.16036 and 0.01066 per site, nucleotide diversity, Pi = 0.00826 and Tajima's D = 20.83312.

Elytral Reticulation and Phylogeny
Parsimony analysis of the combined data set, including all outgroup taxa, resulted in 13 most parsimonious (MP) trees with a length of 400 steps (uninformative characters excluded), and CI and RI of 0.65 and 0.83, respectively. These fundamental trees are combined into a strict consensus (Fig. 5). The subsequent phylogenetic analysis only including the A. bipustulatus complex and A. nevadensis specimens resulted in four MP trees with a tree length of 48 steps, CI and RI of 0.90 and 0.97, respectively. Main differences between these fundamental trees were the position of three haplotypes (BipSwe2, BipFra3 and BipPor2) and the grouping of the two clades consisting of the BipIce1, BipIce2, BipFra12, BipFra13 and the A. nevadensis, BipFra5, BipIra1, and BipMor1 specimens. The topology and subclades of the strict consensus tree from both analyses were identical, which provides strong evidence that no long branch attraction is affecting the results. The Kishino and Hasegawa test of the four MP fundamental trees of the A. bipustulatus complex resulted in one phylogenetic hypothesis with the lowest likelihood value. This tree is hereafter referred to as the ''best'' tree ( Fig. 6).
According to the Akaike Information Criterion in Modeltest, the Translational model (TIM), with a proportion of invariable sites equal to 0.8199 and a gamma distribution shape parameter of 0.7877 is the best evolutionary model, within the maximum likelihood analysis, given the data only including the A. bipustulatus complex and A. nevadensis specimens (-lnL 2954.8911 and AIC score 4125.7822). Nucleotide frequencies; (A) 0.3368, (C) 0.1501, (G) 0.1214 and (T) 0.3917 and evolutionary model parameters; A-C: 1.0000, A-G: 11.6809, A-T: 0.1441, C-G: 0.1441, C-T: 6.3908 and G-T: 1.0000. The maximum likelihood analysis resulted in six trees after 75166 rearrangements tried from a starting Neigbour Joining tree. From these trees the Kishino-Hasegawa (KH) test selected significantly one tree (ln 348.43962) to be the most likely representing the evolutionary history of the data (Fig. 7).
Two main evolutionary lineages are present within the ''best'' MP tree (Fig. 6) and the ML tree (Fig. 7), hereafter referred to as group I and II. These two groups are separated by 10 unambiguous unique nucleotide transformations, which result in both a 100% bootstrap (BS) and Jackknife support (JS). Group I includes specimens from two haplotypes from Russia, two from Sweden and one from France, whereas all other specimens from Sweden, France, Spain, Island, Italy, Morocco, and Iran belong to group II. The ''best'' MP tree (Fig. 6) and the ML tree (Fig. 7) differ chiefly in the position of the BipSwe2 haplotype and the presence or absence of the subclade including the haplotypes BipIce1, BipIce2, BipFra12 and BipFra13, and the subclade including BipEsp4 and BipEsp5. These differences do not change the evolutionary interpretation of the morphological results of this study, however. Both kinds of phylogenetic analyses of the A. bipustulatus complex clearly demonstrate a multiple independent origin of different reticulation patterns. This pattern is evident in the kiesenwetterii form (B type) that is found in several of the monophyletic groups in the MP strict consensus tree together with haplotypes of the bipustulatus form (A type) (Fig. 5 and Fig. 7). The falcozi form also displays a similar pattern and is found within two well-supported subclades (71 and 63% BS, respectively) within the MP analysis, and basally of group II and in one of its subclades in the ML (75% JS). The dolomitanus form (B type) from Italy (BipIta1) forms a strongly supported subclade (96% BS and 99 % JS) together with the bipustulatus form (A type) from France (BipFra6) and Spain (BipEsp3). Moreover, the grouping of bipustulatus, solieri and kiesenwetterii as sister forms within several well supported subclades (BS and JS 53-93%) strengthens the argument of adaptation to altitudinal environments.
Forcing the kiesenwetterii, falcozi or solieri forms on the ''best'' MP fundamental tree (Fig. 6) to be monophyletic resulted in a dramatic increase in three length of 71.1 % compared to the fundamental trees. Still an even larger increase (83.3%) is seen when we add an additional topological constraint assuming A. nevadensis as sistergroup.
Agabus nevadensis is deeply nested in group II, within the A. bipustulatus complex (Fig. 5 and Fig. 7). However, the allozyme differences seen in the population genetic analysis imply that the adaptation process can lead to assortative mating and the evolution of a reproductive barrier.

Discussion
We have studied the phylogeographic history and genetic variation of different reticulation patterns across the A. bipustulatus complex and A. nevadensis distribution in the West Paleartic in order to understand the complex variation in elytral reticulation that occurs in both sexes at high altitudes in Central and South Europe in contrast to the more homogeneous reticulation pattern normally observed at lower altitudes across the total distribution area.
We found that the different secondary reticulation patterns in A. bipustulatus have evolved independently at several times. The phylogenetic pattern shown here within the ''best'' MP fundamental and ML tree describes the A. bipustulatus complex as a deep gene tree with two major lineages broadly sympatric, between specimens with different reticulation patterns (Fig. 6). This reflects a species with a large evolutionary effective population size and high gene flow, with a recent admixture of lineages that have diverged allopatrically [54]. The geographic mix of sampled specimens within the recognised monophyletic group I and II, of A. bipustulatus indicates that the Pyrenees and Alps have not acted as dispersal barriers in the interglacial periods, which is argued to be a major potential isolation barrier for the Iberian water beetles [31]. Instead specimens of A. bipustulatus have been able to disperse more or less randomly into different areas during the climatic oscillations that have occurred during the late Pleistocene [55]. In addition, specimens from both groups I and II co-occur at the Stintbä cken locality in north Sweden. There is therefore no reason to believe that the deep gene tree represents two different monophyletic genetically isolated entities, since the genetic structure of the Stintbä cken population does not differ from other A. bipustulatus populations across its distribution area [56].
Our results also indicate the presence of different selective regimes at different altitudes. The lack of significant genetic differences among A. bipustulatus populations representing combinations of reticulation patterns indicates that the historical legacy of a population does not hinder gene flow [57][58]. However, decreasing mean heterozygosity and increased population substructure (F IS ) at different altitudes regardless of reticulation patterns shows that the selective process, is strong enough to alter/ change the genetic structure of populations ( Fig. 4 and Table 4). This genetic pattern is congruent with earlier findings in north Scandinavia [27,56]. Individual response in secondary reticulation at high altitude is observed as widening of the longitudinally elongated secondary meshes and increased coverage of isomorphic meshes from the apex up to half of the elytral length. Males are more affected than females ( Table 3). The occurrence of especially secondary reticulation type C and G specimens in France and Italy at very low frequencies together with other reticulation patterns outside the Vall de Boí in Catalunya (Table 3) implies that crosses between reticulation patterns may give rise to patterns not observed within parents seen as increased reticulation diversity at higher altitude. The geographical and altitudinal distribution of the reticulation types and the genetic variation within the bipustulatus-complex suggest a complex evolution. The enlarged isomorphic meshes of the secondary reticulation at high altitude in females and the evolution of wider male tarsi and attached suckers in the presence of modified reticulated females could indeed contribute to the arms-race between sexes proposed by Howe [13]. An alternative explanation for the increased variance in female reticulation at high altitude could be that it arises as a side effect of the greater variation of elytral reticulation in males in these populations.
The high altitudinal populations of A. nevadensis stand out in comparison to the genetic variation found among populations of A. bipustulatus and their phylogenetic position deeply nested within the A. bipustulatus group II ( Fig. 6 and 7). The large similarity of the primary and secondary reticulation in A. nevadensis to that of the kiesenwetterii form (type B, Table 2) suggests that the initial altitudinal adaptation follows the same selection process as the other reticulation patterns in A. bipustulatus (type A). However, the significant differences in allele frequencies between A. nevadensis in Sierra Nevada and A. bipustulatus populations in Spain, France, Italy and Sweden at different altitudes strongly indicate that A. nevadensis has relatively recently became reproductively isolated in situ [59][60][61].
Drotz et al. [27] documented that selection acted on a-Glycerophosphate dehydrogenase (a-Gpdh) in A. bipustulatus between the valley floor and above the tree line in Scandinavia. This selective force is, however, not responsible for the observed reticulation variation seen in both A. nevadensis and A. bipustulatus since it occurs across the species total distribution area and affect all analysed forms of both species regarding their reticulation pattern [56]. Studies of the leaf beetle Chrysomela lapponica could, however, give some clue to a possible selective candidate. This species displays a patchy distribution of elytral colour variation with large dark black spots in northern Europe and high mountainous regions in Western Europe. In central Europe, brightly red coloured beetles are more common at low altitude. Gross et al. [62] demonstrated empirically that the increase in melanic elytra in northern Europe and high altitude represent a selective thermal advantage at low temperatures in relation to the brightly coloured central European form. Here the size of the dark black spots regulates the body temperature. Similar demographic distribution and colour morph variation is seen in Coelophora inaequalis, Harmonia axyridis and Oreina sulcata [63][64][65]. The colour variation within these species is assumed to be driven by solar UVb radiation and interestingly solar UV-b radiation normally increases with increasing altitude. This implies that the strongest selective effect within a geographical area is overlapping with our observed deviation in reticulation pattern of A. bipustulatus and A. nevadensis. The Alps, Pyrenees and Sierra Nevada display the highest UV levels in Europe as a consequence of their high altitude, snow-covered surface and low aerosol levels [66]. In addition, the effect of UV radiation in aquatic environments is largest in shallow mountain lakes with high incident flux and deep penetration of UV radiation [67]. This habitat is often associated with A. bipustulatus species at high altitude [68]. Here different reticulation patterns may reflect UV radiation differently where more impressed primary reticulation can lead to increased body temperature [69].
Further research is needed in order to fully understand the interplay between the selective forces that clearly plays a part in forming the reticulation pattern observed within the A. bipustulatus complex and their combined importance for adaptive/sympatric speciation. One way to do this could be to address how the relative strength of sexual and natural selection changes at different altitudes, by analysing male reticulation and suction cup variation in relation to the female reticulation patterns at different altitudes. This, however, is beyond the scope of this paper.