Erebia euryale finds its ecologic optimum in edges and clearings of the timberline forest, irrespective of the forest association (conifer or beech forest), or the underlying bedrock (calcareous or siliceous). Above the timberline, it can be abundant in the green alder shrubs (Alnetum viridis) and in the alpenrose heath vegetation (Rhododendro–Vaccinietum). Locally, it descends into the montane zone or ascends into the alpine meadows (Gonseth 1987; Sonderegger 2005). The larvae feed on grasses, their development is biennial (Sonderegger 2005; Klecková et al. 2015).

Adult individuals were collected at 16 localities (Fig. 1 and Table 1). As an outgroup, two males were sampled of E. ligea (Linnaeus, 1758), the nearest relative of E. euryale (Dincă et al. 2010; Peña et al. 2015; Dincă et al. 2021). Netted individuals were paralysed with ethyl acetate vapour, and two legs were transferred to ethanol 96% before death. Specimens included in this study are all registered on the Barcode of Life Database (Ratnasingham and Hebert 2007) with detailed information on collecting localities and identification (dataset https://doi.org/10.5883/DS-EURYALE) Vouchers are deposited in the collection of Naturalis Biodiversity Center (RMNH; Leiden, the Netherlands). Additional COI sequences were retrieved from BOLD and GenBank, from localities not represented in our samples. These covered bp 1–658 (Dincă et al. 2010), bp 668–1483 (Vila et al. 2011) or bp 305–1475 (Paučulová et al. 2017) (Table 1, Suppl. material 1: Table S1).

Genomic DNA was extracted from the legs with a Macherey–Nagel NucleoMag 96 Tissue magnetic bead kit on a Thermo Fisher KingFisher flex system. Polymerase chain reaction (PCR) was used to amplify a fragment of 1,496 base pairs (bp) of the Cytochrome C Oxidase subunit I gene (COI) from the mitochondrial genome, in two parts. The first part of 658 bp, near the 5’P side of the gene, also known as the DNA barcode region (Hebert et al. 2003) was amplified using a 1:1 primer mix of LepF1 (5’ ATT CAA CCA ATC ATA AAG ATA TTG G 3’) and LCO1490 (5’ GGT CAA CAA ATC ATA AAG ATA TTG G 3’) as forward primers, and a 1:1 primer mix of LepR1 (5’ TAA ACT TCT GGA TGT CCA AAA AAT CA 3’) and HCO2198 (5’ TAA ACT TCA GGG TGA CCA AAA AAT CA 3’) as reverse primers (Folmer et al. 1994; Doorenweerd et al. 2014). The second part that we amplified, from here on indicated as the 3’P fragment, has a 32 bp overlap with the 5’P fragment, and a target length of 870 bp. This fragment was amplified with forward primer SeqIntCOIf (5’CWT CWT TTT TTG AYC CAG CWG GAG 3’) and reverse primer LepLEUr (5’ CCA TTA CWT ATA RTC TGC CAT ATT 3’) (Vila and Björklund 2004). M13 forward (5’ TGT AAA ACG ACG GCC AGT 3’) or reverse (5’ CAG GAA ACA GCT ATG AC 3’) tails were attached to the forward and reverse primers respectively, and used for sequencing. Each PCR reaction included 40 cycles, consisting of 3 minutes initial denaturation at 94 °C, 15 seconds cycle denaturation at 94 °C, 30 seconds cycle annealing at 50 °C, and 40 seconds cycle extension at 72 °C. A final extension at 72 °C for 5 minutes completed the reactions. Bidirectional Sanger sequencing was outsourced to BaseClear, Leiden, The Netherlands. The resulting chromatograms were checked for quality and congruence in Geneious R6.1.8 and the resulting sequences were managed using VoSeq 1.7.4 (Peña and Malm 2012) prior to uploading to BOLD (www.barcodingoflife.org) and Genbank (www.ncbi.nlm.gov/Genbank) (accessions MT762427–MT762602).

Specimens with missing data were removed, leaving 158 E. euryale samples for haplotype analyses. Median-joining haplotype networks were constructed using NETWORK v10.2.0.0 (www.fluxus-engineering.com), based on different sequences: (i) Our main dataset of 1495 bp, i.e., the combined 3’P and 5’P fragments. (ii) Two subsets, of the 3’P and 5’P fragments separately, both expanded with sequences mined from GenBank. (iii) A set of sequences covering bp 305–1475, retrieved from GenBank. Because the mined 3’P sequences were 816 bp long instead of 838, we truncated our 3’P sequences to bp 816 bp. The excised nucleotides contained no mutations. Because the 5’P fragment in ssp. syrmia was invariable (N=17) and obviously did not contribute to the variation, we concatenated the mined 3’P and 5’P sequences of this subspecies and included them into our 3’P+5’P dataset, which we accordingly truncated to 1,474 bp.

From the main dataset we inferred a maximum-likelihood tree using IQ-TREE v1.6.10 (Nguyen et al. 2015) with 5,000 approximate likelihood ratio test repetitions (flag -alrt) and 5,000 rapid bootstraps (flag -bb). The IQ-Tree integrated ModelFinder selected TN+I+F as the best fit model according to the Bayesian information criterion, but all models scored very similar. For Bayesian tree inference and molecular dating we used RevBayes v1.0.12 (Höhna et al. 2016). We calibrated the root of the tree, the split between E. euryale and E. ligea, at 3.49 My to match the estimate that was based on a four gene Erebia wide phylogeny (Peña et al. 2015). A study that used a supertree approach with largely the same genetic data for Erebia but different calibration points, estimated the split between E. ligea and E. euryale very similarly, at 3.58 My (Wiemers et al. 2020). We ran two independent MCMC chains in parallel for 20,000 generations, both with the GTR+G substitution model, and used Tracer v1.7.1 (Rambaut et al. 2018) to verify that this resulted in sufficiently large estimated sampling sizes (>>200).

Reconstructing phylogeography requires insight in past area shifts. Because Würm glacial refugia form the bridge between present and pre-Würm distribution patterns, we primarily examined how the present distributions line up with potential Würm refugia. Six regions are broadly recognised as potential Würm glacial refugia: (i) In the westernmost (French pre-Alps, south of the Isère glacier) and in (ii) in the easternmost parts of the Alps (south of the Enns glacier), extended areas of unglaciated lower mountains acted as large scale refugia to alpine organisms (Schönswetter et al. 2003; Buoncristiani and Campy 2004; van Husen 2004; Schönswetter et al. 2005). (iii) Along the southern border of the Alpine arch, habitable areas were present on the south exposed slopes. These were unglaciated (Penck and Brückner 1909; Geologische Bundesanstalt (AT) 2013), and foothill topography ensured wind shelter, increased insolation, and moist conditions by orogenic precipitation and summery melt water supply, resulting in tempered micro habitats (Avigliano et al. 2002; Kropf et al. 2002; Ravazzi et al. 2004; Monegato et al. 2007). A key feature of this habitable belt is that it was heavily compartmentalized by north–south running glacial valleys (Penck and Brückner 1909; Castiglioni 2004). These valleys acted as habitat interruptions, impeding lateral dispersal. (iv) Along the northern Alpine border, a lobed ice front composed of piedmont glaciers in lateral contact, shovelled far into the foreland (Stojakowits et al. 2020). Some persisting incisions in this ice front (Keller and Krayss 2010; Reitner 2011; Heß 2013; Graf et al. 2015) offered suitable refugial conditions to sub alpine organisms (Schönswetter et al. 2005). (v) Adjacent to the eastern Alps, the Pannonian Basin has been evidenced to have had taiga like forests and even patches of temperate deciduous woodland during the Würm glacial (Rudner and Sümegi 2001; Willis and van Andel 2004; Feurdan 2005; Janovska and Pokorný 2008; Kuneš et al. 2008; Juřičková et al. 2014). (vi) The southern Pyrenean foreland offered favourable conditions even during pleniglacial periods (González-Sampériz et al. 2005). This goes to a lesser degree for the western and eastern edges of the Pyrenees (Calvet 2004). The northern Pyrenees were far less hospitable (Heinz and Barbaza 1998), and treeless (de Beaulieu et al. 1994).

For the nomenclature of Pleistocene warm and cold phases we follow the Marine Isotope Stages (MIS) system (Martinson et al. 1987; Lisiecki and Raymo 2005) (Appendix 1: Fig. A1). The main Würm glacier advance (the Last Glacial Maximum, LGM) corresponds to MIS2, the main Riss advance to MIS6. We avoid the terms Mindel and Günz, because there is no consensus on their correlation with MIS stages (cf. Head and Gibbard 2005; Häuselmann et al. 2007). Both probably cover more than one climatic cycle (Kukla 2005). Nomenclature of mountain stocks is according to Marazzi (2005).

The main network, based on 1,474 bp of COI (Fig. 2A) contains fifty-three haplotypes, determined by 51 variable sites (two variants per site), accounting for 3.4% of the sequence. Thirty-two sites were parsimony informative, 13 of which were located on the 5’P fragment and 19 on the 3’P fragment. Twenty-nine substitutions (57%) were transitions. Fifty-one haplotypes (98%) were private to a single subspecies, forty-six (88%) to a single population. Haplotypes are numbered 1 through 53 in Fig. 2A. They appeared grouped in five clusters, designated A through E. Additional haplotypes in the mined sequences of the 3’P fragment are numbered 54 through 59, (inserts in Fig. 2A). The maximum distance (outgroup excluded) was 9 mutational steps (0.61% pairwise sequence difference). The outgroup Erebia ligea (haplotypes 60 and 61) was connected by 20 mutational steps, a pairwise sequence difference of 1.35%. The only loop, regarding haplotype 37, was resolved by breaking the connection to cluster C according to the criteria by Crandall and Templeton (1993). The 5’P fragment (Fig. 2B) was hardly discriminating. In the network based on bp 305–1,475 (Fig. 2C), the Sudeten and Carpathian samples retrieved from GenBank appear all in cluster E.

Figure 2. 

Median–joining haplotype networks, based on different COI segments. Circle diameters are proportional to the frequency of the haplotype. Red diamonds represent hypothetical haplotypes. A. Based on 168 3’P+5’P sequences of Erebia euryale. Inserts: cluster A and sub cluster E2, based on the 3’P section of the same 168 sequences, together with 77 sequences retrieved from GenBank (Vila et al. 2011). Haplotypes are numbered 1 through 61; B. Based on the 5’P section of the same 168 sequences plus 23 sequences retrieved from GenBank (Dincă et al. 2010); C. Based on bp. 305–1,475 of the same 168 sequences plus 36 sequences mined from GenBank (Paučulová et al. 2017). In Fig. 2C only the part corresponding to cluster E in Fig. 2A is depicted. Colour legends in frame B apply to Figures 2 to 4.

The relationship between subspecies, clusters and haplotypes is shown in Table 2. The clusters in the network turned out largely congruent with morphologic subspecies. Cluster A (Fig. 2A) includes all haplotypes of the western extra-Alpine populations: the Jurassian ssp. tramelana, the Pyrenean ssp. pyraenaeicola and the Cantabrian ssp. cantabricola.

Subspecies adyte makes up cluster B, but several adyte specimens are placed in cluster E.

Cluster C contains all individuals of ssp. pseudoadyte. It has the highest genetic diversity of all subspecies with 11 haplotypes found in two localities.

Cluster D contains only haplotypes of ssp. kunzi, but some kunzi are assigned to cluster E. Western and eastern kunzi are not separated. The Monte Cavallo population (sample 11) builds a separate sub cluster D1, four mutational steps apart from D2 (Fig. 2A), but morphologically indistinguishable from it.

Cluster E is private to four subspecies (euryale, syrmia, isarica and ocellaris), and contains individuals of two subspecies that have their own private clusters (adyte and kunzi). In the first group there is a clear structuring: ssp. syrmia is found around the central axis of the network, from which four subclusters, E1 through E4, branch off (Fig. 2A, C). The sub-clusters E2 and E3 are private to ssp. isarica, E4 to ssp. euryale. Ssp. ocellaris occupies E1, but is also represented in E3, together with isarica, with which it has a documented hybrid zone (Cupedo 2010). The proportion of individuals with the shared haplotype increases with decreasing distance to the hybrid zone (Table 2). The genetic dichotomy in ssp. isarica (E2 and E3) coincides with a geographic-morphologic sub-division. E3 contains the western and northern Alpine populations, E2 the eastern Alpine populations. These groups differ in wing shape and wing pattern (pers. obs. FC).

As to the subspecies with a different private cluster, we suspect that their cluster E haplotypes are due to introgression or incomplete lineage sorting. In sample 4 there is 100% mismatch of phenotype and haplotype. Phenotypically it is an adyte enclave in isarica territory (Cupedo 2010). The haplotype matches the surrounding isarica populations. This fits with previous studies in which gene-flow between these subspecies was not reflected in morphology (Geiger and Rezbanyai 1982). In ssp. kunzi, on the other hand, the cluster E individuals were found in the most peripheral samples (7 and 10, Fig. 1), far from the boundary with ocellaris, with which they share haplotype 51 (Fig. 2A). This argues against a recent introgression. The fact that both western and eastern kunzi are involved, suggests an introgression event predating their separation.

Bayesian Inference (BI) and Maximum Likelihood (ML) analyses based on the 1,474 bp alignment (Fig. 3) show similar topologies for the haplogroups A through D, but the position and relationships of haplogroup E are not well supported in any of the trees. The clade comprising ssp. adyte and the western extra Alpine subspecies (the network clusters A+B), is well-supported in both trees. The western extra Alpine subspecies are jointly monophyletic with good support, a monophyletic ssp. adyte is not supported. There is medium support for a monophyletic ssp. pseudoadyte, private cluster C in the network. The clades D1 and D2 are well supported in both trees, the monophyly of D1+D2, however, is not well supported in the BI analysis and only moderately supported by ML. The monophyly of individual clades E1 and E2 is also well supported in both analyses, E3 somewhat less so in BI, but their positions in the trees are not well resolved. Overall, the BI with molecular dating suggests that the diversification into the present subspecies took place between ~0.75 and ~2 My ago.

Figure 3. 

Phylogenetic trees based on the haplotypes in Fig. 2A, with clusters and sub-clusters collapsed into triangles. Colour codes as in Fig. 2, cluster labels as in Fig. 2A. A. Bayesian Inference tree. Support values on the branches indicate the Bayesian posterior probabilities. Blue node bars indicate the 95% posterior density for the Bayesian age estimates. Insert: estimates of the divergence time from the RevBayes approach, based on a calibrated split at 3.5 My BP between E. euryale and E. ligea (Peña et al. 2015). The time scale is in My ago. B. Maximum Likelihood tree. Support values on the branches indicate bootstrap support (5,000 replicates).

The actual distribution pattern results from Holocene range extensions, originating in the Würm glacial refugia. Overlaying the potential Würm refugia with the current subspecies distributions shows a clear matching pattern.

Ssp. adyte occupies the entire southwestern Alps (Fig. 1). Its distribution area stretches between two known biogeographical borders (Cupedo and Doorenweerd 2020): the Isère valley in the west and the Valtellina (the Adda valley) in the east, suggesting survival in refugium (i) and in the western compartments of refugium (iii).

The territory of ssp. pseudoadyte is sharply delimited by the valleys of Adda and Adige (Fig. 4), matching a single compartment of refugium (iii). Its postglacial northward expansion was halted at the Adda and Adige valleys, which currently delimit its distribution.

Figure 4. 

The southern Alps east of Lake Como, showing the relationship between the current distribution (coloured areas) and hypothesized refugial areas (coloured dots) of the subspecies of E. euryale during the maximum Würm glaciation (MIS2). Dotted area – transitional populations isarica/ocellaris. Colour legends as in Fig. 2.

As to ssp. kunzi, all former peripheric nunataks between Val d’Adige and Valcellina, belonging to refugium (iii), harbour relict populations of this ssp. (Cupedo 2010; Bonato et al. 2014): the Monti Lessini, Monte Grappa, Col Visentin and Monte Cavallo (Fig. 4). Sub-cluster D1 to date is confined to its refugial massif, the Monte Cavallo, and is separated from D2 by the Piave valley (Fig. 4). Only sub-cluster D2, living across the Piave valley, could invade the southern Dolomites, where it failed to cross the deeply incised Cordevole valley (Fig. 4). Eastern and western D2, only differing in wing pattern, are separated by the Cismon-Brenta valley.

The ssp. ocellaris territory is considerably larger than that of ssp kunzi, but it includes only a short section of the southern pre-Alpine chain (Figs 1, 4). Its refugium along the Alpine border most probably stretched from the Valcellina to the Tagliamento valley, matching the eastern part of refugium (iii).

Ssp. isarica has by far the largest distribution of the Alpine subspecies (Fig. 1). We suppose that the two COI sub clusters represent populations originating from different glacial refugia. The current territory of cluster E2 comprises the vast refugial area (ii) on the eastern edge of the Alps (Fig. 4), whereas cluster E3 suggests refugia along the northern outskirt of the Alps (refugium iv), the exact locations of which remain unknown.

The five-fold cluster structure of the network is superimposed by a geographic tripartition. We distinguish (i) a southern Alpine group, consisting of the clusters B, C and D, (ii) a circum–Pannonian group, coinciding with cluster E, in which Carpathian populations are combined with northern and eastern Alpine ones and (iii) a western extra Alpine group, equalling cluster A.

Ad (i). In the southern Alpine group the subspecies adyte, pseudoadyte and kunzi occupy adjacent territories. They share a sharply defined, coherent territory and a considerable genetic variation. The pairwise sequence difference within the clusters ranges from 0.68% to 1.15%. This suggests a long-term isolation within the actual distributional borders. Most probably the increasing incision of the dividing valleys, by repeated glacial erosion, resulted in a permanent vicariance, enabling their differentiation.

Ad (ii). In the circum–Pannonian group, the core area of cluster E in the haplotype network is occupied by the Southern and Eastern Carpathian ssp. syrmia (Figs 1A, 2A). Sub-clusters radiate from it into the Western Carpathians (ssp. euryale, Fig. 2C) and into the Alps (ssp. isarica and ssp. ocellaris, Fig. 2A).

This division into a southern Alpine group and a circum–Pannonian group suggests an ancient Alpine–Carpathian disjunction. This is in agreement with earlier allozyme data: Schmitt and Haubrich (2008) found the split between the Carpathian and Alpine lineages (in their study: ssp. syrmia and ssp. adyte) the most ancient disjunction in the studied E. euryale populations. Such an Alpine–Carpathian disjunction is typically found in species of humid grasslands of the montane-subalpine levels (Varga and Schmitt 2008), and is associated with a Pannonian refugium (refugium v). The common ancestor could, in theory, have been resident either in the Alps or in the Carpathians. The close genetic affinity between the outgroup E. ligea and the Carpathian samples in the network (Fig. 2A) suggests a Carpathian origin, but the phylogenetic tree lacks the statistical support to corroborate this scenario.

Subsequently, i.e., during subsequent cold periods, the southern Alpine group differentiated into its three lineages (adyte, pseudoadyte and kunzi), and the circum-Pannonian group split up into a Carpathian and an eastern Alpine lineage. This required a second, more recent glacial retreat in the Pannonian plain, again followed by interglacial invasion of both the Alps and the Carpathians. Finally, the Sudeten and Carpathian populations differentiated into what is now known as ssp. euryale and ssp. syrmia, while the Alpine populations differentiated into today’s ssp. ocellaris, eastern isarica and western isarica. We presume that this occurred during MIS2, when the refugia of the Alpine populations were localised in the Alpine periphery rather than in the Pannonian plain (Fig. 4), because the Alpine glacier extension was far less than during the foregoing glaciations (Penck and Brückner 1909; Castiglioni 2004; van Husen 2004).

Ad (iii) The western extra-Alpine group. Both the network and the ML tree suggest that it derived from cluster B (Figs 2A, 3B). Its star-like configuration, with many low frequency haplotypes derived from a most common one, is typical of a relatively recent demographic expansion after a genetic bottleneck (Rogers and Harpending 1992; Avise 2000; Hwang and Cho 2018). This bottleneck probably was associated with its separation from the Alpine cluster B (Fig. 2A). Because an Alpine–Pyrenean disjunction required sub alpine conditions in the in-between Rhone basin, it must have been associated with a glacial period. Because of the identical genetic content of the Pyrenean and Jurassian populations, we exclude two independent colonisation events. We assume that either the Pyrenees and the Jura were colonised simultaneously, or the Jura was colonised secondarily from the Pyrenees.

However, nuclear data suggest a somewhat different scenario. Schmitt and Haubrich (2008) found the Pyrenean population genetically closer to isarica than to adyte. Although their data refer only to eastern isarica, this suggests a connection to the northern, rather than to the southern Alps. This mismatch between mitochondrial and nuclear data calls for further study on the origin of the western extra-Alpine group of E. euryale. In general, mitochondrial and nuclear DNA need not have responded similarly to past biogeographic events, and their patterns of variation provide complementary data on the phylogeographic history (Hinojosa et al. 2019). We expect additional genomic data to enhance the insight in the origin of the western extra-Alpine populations.

Some morphologic taxonomic characters have been shown phylogenetically relevant, provided they are combined with molecular data. (Wahlberg and Nylin 2005; Simonsen et al. 2006; Pisani et al. 2007; Stephanović et al. 2016). In our study, the three morphologic subspecies that make up the southern Alpine group (adyte, pseudoadyte and kunzi) appeared congruent with well separated haplotype clusters (Fig. 2A). These three subspecies (together with the Apennine ssp brutiorum which is not included in this study), share a distinctive character: they have white pupiled apical ocelli on the forewing upper side. All remaining subspecies (including the boreal/tundral populations, not studied either), have uniform black ocelli. The pupiled ocelli appear to be a synapomorphy of the southern Alpine group. This supports the conclusion that the separation of the southern Alpine group was a basal split in the history of Erebia euryale. It further suggests that the Apennine populations are derived from the southern Alpine group, and that the boreal/tundral populations are derived from the circum-Pannonian group. Both hypotheses are consistent with geography. Moreover, it supports the alternative, allozyme based scenario for the origin of the western extra-Alpine group, in which the apical ocelli are black.

The Bayesian estimate of divergence times, using an external calibration, dates the differentiation of the crown group roughly 0.75 My to 2 My ago, i.e., in the early Pleistocene (Appendix 1: Fig. A1). However, the correlation of lineage splits to glacial conditions provides an alternative calibration for the age of splitting events (Schmitt et al. 2016; Hinojosa et al. 2018). For the relatively young Pyrenean disjunction we assume at least a MIS6 age, because there is some evidence for more recent (i.e., MIS2) differentiations within this cluster: (i) genetic differences between populations at both sides of the Cantabrian watershed are a possible indication of two Würm glacial refugia in the Cantabrian region (Vila et al. 2011). (ii) Likewise, the morphological differentiation of the Pyrenean populations into ssp. antevortes in the north-western Pyrenees (not investigated here) and ssp. pyraenaeicola in its southern and eastern parts, must be younger than the Alpine–Pyrenean disjunction.

Because the differentiation of the clusters B, C and D in their southern Alpine refugia predates the Alpine–Pyrenean disjunction (Fig. 3), it is at the latest of MIS10 age. The initial Alpine–Carpathian disjunction, necessarily associated with a severe cold phase, thus most probably is not younger than MIS12. We emphasise that the ages of disjunctions reconstructed this way are always minimal ones and splitting events might prove considerably older.

The dependency of major cold phases, however, sets a maximum age too. It has been shown that tributary glaciers in the eastern Alps did not reach the main valleys during minor cold periods like MIS4 (van Andel and Tzedakis 1996; van Husen 2004), therefore subalpine organisms were not dispelled into refugial isolation. This was only the case during the more severe cold phases. The very first of those major ice ages was MIS22, ca 870–880 ka ago (Ehlers and Gibbard 2007) (Appendix 1: Fig. A1). All foregoing cold phases were in the order of MIS4 or less. Consequently, it is unlikely that the first Alpine–Carpathian disjunction took place more than 1 My ago. This conflicts with our Bayesian dating (Fig. 3A), secondary calibrated with Peña et al. (2015). It is in line, though, with recent molecular studies on the E. tyndarus group (Schmitt et al. 2016) and the E. epiphron–orientalis group (Hinojosa et al. 2018). Both calculated divergence times, in the Pliocene–Pleistocene period, that are considerably younger than those inferred from secondary calibrations based on the external calibration from Peña et al. (2015).