# Cross-continental phylogeography of two Holarctic Nymphalid butterflies, Boloria eunomia and Boloria selene

Jana Maresova, Jan Christian Habel, Gabriel Nève, Marcin Sielezniew, Alena

Bartonova, Agata Kostro-Ambroziak, Zdenek Faltynek Fric

► To cite this version:

Jana Maresova, Jan Christian Habel, Gabriel Nève, Marcin Sielezniew, Alena Bartonova, et al.. Cross-continental phylogeography of two Holarctic Nymphalid butterflies, Boloria eunomia and Boloria selene. PLoS ONE, 2019, 14 (3), pp.e0214483. 10.1371/journal.pone.0214483 . hal-02086734

HAL Id: hal-02086734

https://amu.hal.science/hal-02086734

Submitted on 21 May 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Cross-continental phylogeography of two Holarctic Nymphalid butterflies, *Boloria eunomia* and *Boloria selene*

Jana Maresova¹,²*, Jan Christian Habel³, Gabriel Neve¹,⁴, Marcin Sielezniew¹,⁵,
Alena Bartonova¹,², Agata Kostro-Ambroziak⁵, Zdenek Faltynek Fric¹

1 Biology Centre of the Czech Academy of Sciences, Institute of Entomology, Ceske Budejovice, Czech Republic, 2 University of South Bohemia, Faculty of Sciences, Ceske Budejovice, Czech Republic,

3 Terrestrial Ecology Research Group, Department of Ecology and Ecosystem Management, School of Life Sciences Weihenstephan, Technische Universität München, Freising, Germany, 4 IMBE, Aix Marseille Université, CNRS, IRD, Avignon Université, Marseille, France, 5 Laboratory of Insect Evolutionary Biology and Ecology, Institute of Biology, University of Białystok, Białystok, Poland

* maresovajana2@gmail.com

<figure><img src="image_2.png" /></figure>

## Abstract

Pleistocene glaciations had significant effects on the distribution and evolution of species inhabiting the Holarctic region. Phylogeographic studies concerning the entire region are still rare. Here, we compared global phylogeographic patterns of one boreo-montane and one boreo-temperate butterflies with largely overlapping distribution ranges across the Northern Hemisphere, but with different levels of range fragmentation and food specialization. We reconstructed the global phylogeographic history of the boreo-montane specialist *Boloria eunomia* (n = 223) and of the boreo-temperate generalist *Boloria selene* (n = 106) based on mitochondrial and nuclear DNA markers, and with species distribution modelling (SDM). According to the genetic structures obtained, both species show a Siberian origin and considerable split among populations from Nearctic and Palaearctic regions. According to SDMs and molecular data, both butterflies could inhabit vast areas during the moderate glacials. In the case of *B. selene*, high haplotype diversity and low geographic structure suggest long-lasting interconnected gene flow among populations. A stronger geographic structuring between populations was identified in the specialist *B. eunomia*, presumably due to the less widespread, heterogeneously distributed food resources, associated with cooler and more humid climatic conditions. Populations of both species show opposite patterns across major parts of North America and in the case of *B. eunomia* also across Asia. Our data underline the relevance to cover entire distribution ranges to reconstruct the correct phylogeographic history of species.

<figure><img src="image_3.png" /></figure>

### OPEN ACCESS

**Citation:** Maresova J, Habel JC, Neve G, Sielezniew M, Bartonova A, Kostro-Ambroziak A, et al. (2019) Cross-continental phylogeography of two Holarctic Nymphalid butterflies, *Boloria eunomia* and *Boloria selene*. PLoS ONE 14(3): e0214483. https://doi.org/10.1371/journal.pone.0214483

**Editor:** Łukasz Kajtoch, Institute of Systematics and Evolution of Animals Polish Academy of Sciences, POLAND

**Received:** December 4, 2018

**Accepted:** March 13, 2019

**Published:** March 26, 2019

**Copyright:** © 2019 Maresova et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All sequence files are available from the Genbank database (accession number(s) are provided in S1 Table).

**Funding:** The study was granted by the University of South Bohemia (168/2013/P) (ZFF) and (152/ 2016/P, 038/2019/P) (JM) and the Fulbright Association (ZFF). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

## Introduction

The Quaternary period is characterized by severe fluctuations of long glacial and shorter interglacial periods. These climatic changes caused repeated latitudinal expansions and retractions
**Competing interests:** The authors have declared that no competing interests exist.

of species' distributions [1,2]. According to the "classical" refugia scenario, thermophilic species survived glacial periods in geographically distinct southern refugia, disappeared from major parts of their northerly distribution, and re-expanded again in the wake of postglacial warming [3,4]. Cold-adapted species reacted contrary, with shifts into higher latitudes and altitudes during the warmer interglacial stages [5,6]. However, a number of recent studies showed that the situation is much more complex, emphasizing the importance of the extra-Mediterranean refugia for both temperate and cold-adapted species [7].

Generally, cold-adapted species presently inhabit disjunct areas, consisting of a rather continuous circumpolar distribution and isolated patches restricted to higher elevations in southern mountains [8]. Depending on specific habitat association, arctic-alpine species would inhabit open northern treeless tundras plus mountains above the tree line, while boreo-montane species would occur in northern and mountain taigas and in bogs [8]. For all the northern species, Schmitt [9] proposed two biogeographic scenarios: 1) interglacial population isolation persisted during glacial eras, resulting in deeper genetic splits; or 2) continuous distribution throughout the ice ages with subsequent gene flow, interrupted in the wake of postglacial warming, causing rather weak genetic differentiation.

Although both arctic-alpine and boreo-montane species ranges encompass large parts of the Holarctic region [8], studies on the latter are still underrepresented in the literature on the biogeography of biota during the Pleistocene. Few boreo-montane taxa were studied in the western Palaearctic region [10–15], neglecting the fact that most of them occur across the Northern Hemisphere. Considering the vast extent and largely ice-free status of north-central Asia and north-western North America during the Quaternary, these areas must have played an important role in structuring high-latitude biodiversity and boreo-montane species [16–18]. For instance, the Beringian region might have served as northern refugium for many boreal plant and animal species [19,20]. However, only few studies on the phylogeography of boreo-montane species, including entire distribution ranges, have been published so far (cf. [21–23]).

In this study, we contribute to the knowledge of the phylogeography of two species considering their entire distribution ranges across the Northern Hemisphere. We compare the genetic structures of two Holarctic butterflies, the boreo-montane Bog Fritillary *Boloria eunomia* (Esper, 1800) and the boreo-temperate Small Pearl-bordered Fritillary *Boloria selene* (Denis & Schiffermüller, 1775). Both species show similar distribution ranges, including Asia, Europe and North America [24–26], but the two model taxa differ in respect of their habitat demands and ecological specialisation: *Boloria eunomia* inhabits bogs and fens across the Holarctic region, syntopically with species such as *Colias palaeno* (Linnaeus, 1761) and *Oeneis jutta* (Hübner, [1805]) [27]. *Boloria eunomia* shows a disjunct distribution across central and southern Europe and occur restricted to wet meadows and bogs at higher altitudes, with high abundances of the single larval host plants *Polygonum bistorta* L. and *P. viviparum* L. [28,29]. This butterfly is threatened due to the habitat loss or incorrect habitat management across central Europe [30], but may swiftly enlarge its distribution when the habitats are suitable [15,25]. Towards North and East, *B. eunomia* becomes more broadly distributed and common, as it inhabits extensive habitats, which are both moist tundra and willow seeps as well as raised peat bogs. There, *Viola palustris* L. and *Oxycoccus palustris* Pers. are known as additional food plants [25,31,32]. In contrast, the more generalist *B. selene* uses a broader variety of different habitats such as forest clearings, wet meadows, bogs and marshy areas near lakes. The latter species is distributed across the entire Holarctic region, but in contrast to *B. eunomia*, it does not reach the Siberian arctic and it occurs also at lower altitudes and latitudes [33]. Larval food plants are several *Viola* spp. [34–36].
In our study, we applied molecular analyses (mitochondrial and nuclear genes) and performed species distribution models (SDMs) to reconstruct the phylogeography of these butterflies. We hypothesize that *B. eunomia*, inhabiting bogs and humid meadows, represent a strong genetic differentiation due to its rather patchy occurrence during glacial and interglacial stages. On the other hand, we assume that the *B. selene* populations are not or only marginally genetically differentiated among each other as this species probably occupied the same open habitats during both eras, and thus population disjunction occurred only recently.

# Material and methods

## Sampling

We sampled 124 individuals of *Boloria eunomia* and 60 individuals of *B. selene*, covering the whole Holarctic region. The samples were obtained by the authors and from private and public museum collections. A further 99 sequences of *B. eunomia* and 46 of *B. selene* were added to the dataset from GenBank (http://www.ncbi.nlm.nih.gov/) (S1 Table, including specimen ID numbers, GenBank accession numbers, information about specimen repository).

## Acquisition of genetic data

DNA was extracted from legs or the abdomen using the Genomic DNA Mini Kit—Tissue (Geneaid) following the manufacturer's protocols. One fast mitochondrial (*cytochrome c oxidase subunit I*, COI) and two more conservative nuclear (*arginine kinase*, ArgKin; *wingless*, WG) genes were amplified by Polymerase Chain Reactions (PCR) in 20 µl volume (10 µl PPP Mastermix Top-Bio, 6.9 µl PCR H2O, 0.8 + 0.8 µl primers, 1.5 µl DNA). Two forward-reverse primer pairs were used for COI: LCO/HCO and Ron/HCO; for ArgKin, one pair: ArginineF/ArginineR; and for *wingless*, one pair: LepWG1/LepWG2 [37,38]. The universal tails T7Promoter and T3 were attached to all primers. The thermal cycling profile was 95 °C for 5 min, 94 °C for 30 s, 50 °C for 30 s, 72 °C for 90 s, for 36 cycles, and final extension 72 °C for 10 min. Sequencing was provided by Macrogen Inc. on ABI3730XL DNA analysers. Sequences were checked manually and aligned in Geneious v. 8.0.5 [39].

## Genetic diversity and population structure

Maximum parsimony haplotype networks were reconstructed to illustrate relationships among haplotypes of *B. eunomia* and *B. selene* using TCS algorithm [40] in PopArt [41]. The population structure was evaluated by spatial analysis of molecular variance (SAMOVA) implemented in SPADS 1.0 [42], with 10,000 simulated annealing steps and 10 repetitions. This analysis allows defining groups of populations without a priori information on population structure. The method uses a simulated annealing algorithm that maximizes the proportion of genetic variance ($\Phi_{CT}$) partitioned among geographically adjacent populations. SAMOVA was performed for each possible number of population clusters (K) from 2 to 20, with $\Phi_{CT}$ calculated for each value of K. To calculate descriptive statistics of population groups defined by SAMOVA, we used Arlequin 3.5 [43,44].

## Tree construction and time estimates

The best substitution models were selected using PartitionFinder 1.1.1 (see Table 1) [45]. Phylogeny and time diversification of *Boloria eunomia* (COI, ArgKin) and *B. selene* (COI) haplotypes were estimated in BEAST 1.8.4 [46]. We used an uncorrelated lognormal relaxed clock with a coalescent constant size tree and normal distribution priors. For model calibration, we used reference points from the family Nymphalidae [47,48] and the genus *Boloria* [49]; fossils
<table><caption>Table 1. Substitution models for partitioned analyses identified by the PartitionFinder.</caption><thead><tr><th>Partitions</th><th></th><th>PartitionFinder</th><th>BEAST</th></tr></thead><tbody><tr><td colspan="4"><em>Boloria eunomia</em> mtDNA (COI)</td></tr><tr><td>p1</td><td>COI_pos1, COI_pos3</td><td>K81UF</td><td>GTR</td></tr><tr><td>p2</td><td>COI_pos2</td><td>TIM + G</td><td>GTR + G</td></tr><tr><td colspan="4"><em>Boloria eunomia</em> nDNA (ArgKin)</td></tr><tr><td>p1</td><td>ArgKin_pos1, ArgKin_pos2</td><td>F81</td><td>HKY</td></tr><tr><td>p2</td><td>ArgKin_pos3</td><td>TVM + G</td><td>GTR + G</td></tr><tr><td colspan="4"><em>Boloria selene</em> mtDNA (COI)</td></tr><tr><td>p1</td><td>COI_pos1, COI_pos2</td><td>HKY + I</td><td>HKY + I</td></tr><tr><td>p2</td><td>COI_pos3</td><td>GTR</td><td>GTR</td></tr></tbody></table>

If the model was not implemented in BEAST, the first available model with the lowest Bayesian Information Criterion (BIC) scores was used.

https://doi.org/10.1371/journal.pone.0214483.t001

of *Boloria* are unknown. We used three calibration points in *B. eunomia*: first, 12.8 ± 2 Mya for the split between outgroups (*B. aquilonaris* (Stichel, 1908), *B. napaea* (Hoffmannsegg, 1804), *B. pales* (Denis & Schiffermüller, 1775)) + *B. eunomia* and a common ancestor of their sister species; second, the split between outgroups and *B. eunomia* 10.2 ± 2.0 Mya; third, the split between populations from Eurasia and North America 2.0 ± 0.5 Mya. In *B. selene*, we used also three secondary calibration points: first, 9.2 ± 2.0 Mya was the split between outgroups *B. improba* (Butler, 1877) and *B. thore* (Hübner, [1803]) + *B. selene*; second, the split between *B. selene* and *B. thore* 5.0 ± 1.2 Mya; third, the split between populations from Eurasia and North America 2.8 ± 0.8 Mya. We specified MCMC chain length of 50,000,000 generations, sampling every 1000 steps. Operators were auto-optimised. The parameter *ucld.mean* was set as lognormal, while other parameters were left in default values. Three independent runs were performed on each dataset. Tracer 1.6 [50] was used to evaluate estimated values and effective sample size (ESS) for each model parameter and to check convergence. All three runs were combined, and a single best tree was identified. The trees were visualized using the software FigTree 1.3.1. We ran the analyses using the Cipres server [51].

## Species distribution model (SDM)

To investigate climatically suitable areas for the species, we build SDM using georeferenced species records for both species from our own or our colleagues' fieldwork, from publicly available data from the Global Biodiversity Information Facility (GBIF; http://www.gbif.org) and from the literature [52–56]. Duplicates were erased. This resulted in a total of 161 records for *B. eunomia* and 280 for *B. selene* covering the whole species ranges. We obtained information on the present climate and for the climate during the Last Glacial Maximum (LGM; 22000 years ago) as the 19 BIOCLIM variables (BIO1–BIO19) described by the community climate system model (CCSM), the Worldclim 1.4 database [57], at 2.5 arc min spatial resolution. We restricted the area to the Northern Hemisphere (35–75°N, –180°W to 180°E). We calculated the SDMs using MaxEnt 3.4.1 [58] via the R environment using packages 'dismo' and 'raster' [59,60]. Since the records tend to be biased towards the accessible locations (cf. [61]), we used bias background correction for road density [Center for International Earth Science Information Network (CIESIN) Columbia University & Information Technology Outreach Services (ITOS) University of Georgia, 2013], which we rasterized in QGIS v.2.14.3 [62]. We tested the records against 10,355 random background points selected by the script 'sampleRast' (available at https://rdrr.io/github/adamlilith/enmSdm/man/sampleRast.html). We obtained the best values of the intrinsic MaxEnt estimates, regularization multiplier (RM) and the response of
feature classes (FC) using the 'ENMeval' R package [63]. We selected the variables with lowest AICc (AIC corrected for small samples) values after K - 1 jackknife data partitioning. We measured the relative importance of the BIOCLIM variables in the MaxEnt using the jacknife test and excluded the variables with negative gain or gain close to zero expressed as the difference between full model training gain and the gain of the model without the variable(s).

## Results

ArgKin was amplified only for *B. eunomia* and WG only for *B. selene*, possibly due to mutations in the primer region. The latter did not show any intraspecific variation throughout the sampled populations (91 specimens: for GenBank accessions numbers see S1 Table). Consequently, subsequent molecular analyses were conducted only with COI for *B. selene* and COI and ArgKin for *B. eunomia*.

### Boloria eunomia cytochrome c oxidase

We obtained 670 bp COI sequences for 223 individuals from 95 localities and identified 46 haplotypes (Fig 1A and 1B). According to the haplotype network (Fig 1B), haplotypes from Canada and USA (H1-7) formed a separate group differentiated from the Palaearctic, characterised by a widespread haplotype (H1) occurring in high frequencies. Within the Palaearctic region, relatively high haplotype diversity was observed, the highest across Russia with affinities to Scandinavia (H8-16, H20-22, H34-46). European populations were considerably differentiated (H17-19, H23-33). The populations from Cantabrian Mts. (H28) and the Pyrenees (H29-30) slightly differed from each other and were simultaneously related to populations from Ardennes and Nièvre (H32-33). Within Central Europe, two lineages were detected in the Alps: H23-24 and H25-27, the latter closely related to the population from the Czech Republic (the Bohemian Forest; Fig 1 and Figure A in S1 Fig). Poland might represent a contact zone for different lineages (H10-16, H18-19 and H39; Fig 1 and Figure B in S1 Fig). Populations from Eastern Europe (H10-14, H16, H18-19), Balkan region (H17) and Ural Mts. (H15) are closely related (Fig 1).

In SAMOVA, the largest decrease of Φ<sub>CT</sub> variance was for *K* = 11 (Φ<sub>CT</sub> = 0.77), distinguishing 11 population groups: (1) A1-CZ: Alps 1 & CZ; (2) A2: Alps 2; (3) CA-PY: Cantabria Mts. & Pyrenees; (4) AR-NI: Ardennes & Nièvre; (5) EE-UR: Eastern Europe & Ural Mts.; (6) CS: Central Scandinavia; (7) NS-SIB1: Northern Scandinavia & Siberia 1; (8) PL: Poland; (9) SIB2: Siberia 2; (10) CAU: Caucasus Mts. and (11) NA: North America (Fig 2A and Figure A in S2 Fig). The total haplotype diversity was high (Hd = 0.92), with the highest values for NS-SIB1 (Hd = 0.87), SIB2 (Hd = 0.83) and CAU (Hd = 0.71). Significant negative values of Tajima's *D* for A1-CZ and NA and Fu's *Fs* for EE-UR, NS-SIB1 and NA indicate recent population expansion, genetic hitchhiking or selection [64,65] (see Table 2).

The first split of the monophyletic *B. eunomia* clade based on BEAST estimate (Fig 1C) occurred between the Nearctic and Palaearctic haplotypes during the Pleistocene (~1.6 Mya), making the North American clade sister to all Eurasian populations. The first split within the Palaearctic region took place ~0.8 Mya, dividing the haplotypes from Central Scandinavia, Russian Far East and Sakhalin regions from the rest. The remaining Eurasian haplotypes further split into two main branches ~0.6 Mya: One branch with haplotypes from Northern Scandinavia, Poland, Omsk, Altai and Sayan Mts., Novosibirsk and Yakutia; the second with haplotypes distributed mostly in Europe (Spain, Andorra, France, Belgium, Alps, Czech Republic, Bulgaria, Poland, Belarus, Ukraine, Baltic States, S Finland, W Russia) and the Caucasus Mts., Ural Mts., Kamchatka and Chukotka. These haplotypes differentiated around ~0.5
<figure><img src="image_4.png" /><figcaption>**Fig 1.** Phylogeographic patterns of *Boloria eunomia* based on the mitochondrial (COI) marker. Haplotypes distribution on the map COI (A). TCS haplotype network illustrating relationships of the 46 COI haplotypes (B); mutations are shown as hatch marks; the colours correspond to the map (A). Chronogram calculated in BEAST based on COI (C); node bars display 95% highest posterior density interval of the molecular clock and branch labels show posterior probabilities.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0214483.g001
<figure><img src="image_5.png" /><figcaption>**Fig 2.** Spatial analyses of molecular variance (SAMOVA) displayed on the maps. Eleven mitochondrial (COI) lineages of *Boloria eunomia* (A); three nuclear (ArgKin) lineages of *Boloria eunomia* (B); six mitochondrial (COI) lineages revealed in *Boloria selene* (C).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0214483.g002
<table><caption>Table 2. Summary of genetic diversity indices and neutrality tests for *Boloria eunomia* (COI, ArgKin) and *Boloria selene* (COI).</caption><thead><tr><th>Group</th><th></th><th>N</th><th>S</th><th>h</th><th>Haplotypes</th><th>Hd ± SD</th><th>π ± SD</th><th>Tajima's D</th><th>Fu's FS</th></tr></thead><tbody><tr><td colspan="10"><em>Boloria eunomia</em> mtDNA (COI)</td></tr><tr><td>1</td><td>Alps 1 & CZ (A1-CZ)</td><td>14</td><td>6</td><td>3</td><td>H25—H27</td><td>0.2747 ± 0.1484</td><td>0.0013 ± 0.0011</td><td>- 1.9589**</td><td>0.5351</td></tr><tr><td>2</td><td>Alps 2 (A2)</td><td>10</td><td>6</td><td>3</td><td>H23, H24</td><td>0.6222 ± 0.1383</td><td>0.0035 ± 0.0023</td><td>0.3710</td><td>0.5351</td></tr><tr><td>3</td><td>Cantabria Mts. & Pyrenees (CA-PY)</td><td>25</td><td>3</td><td>4</td><td>H28—H31</td><td>0.6167 ± 0.0639</td><td>0.0016 ± 0.0012</td><td>0.8924</td><td>0.4115</td></tr><tr><td>4</td><td>Ardennes & Nièvre (AR-NI)</td><td>14</td><td>1</td><td>2</td><td>H32, H33</td><td>0.2637 ± 0.1360</td><td>0.0004 ± 0.0005</td><td>- 0.3414</td><td>0.1857</td></tr><tr><td>5</td><td>Eastern Europe & Ural Mts. (EE-UR)</td><td>50</td><td>10</td><td>9</td><td>H10—H17, H39</td><td>0.6686 ± 0.0706</td><td>0.0018 ± 0.0013</td><td>- 1.2880</td><td>- 3.1749*</td></tr><tr><td>6</td><td>Central Scandinavia (CS)</td><td>10</td><td>1</td><td>2</td><td>H35, H37</td><td>0.2000 ± 0.1541</td><td>0.0003 ± 0.0005</td><td>- 1.1117</td><td>- 0.3393</td></tr><tr><td>7</td><td>N Scandinavia & Siberia 1 (NS-SIB1)</td><td>25</td><td>12</td><td>11</td><td>H8, H9, H34, H38, H40—H46</td><td>0.8733 ± 0.0496</td><td>0.0039 ± 0.0024</td><td>- 0.6109</td><td>- 3.5831*</td></tr><tr><td>8</td><td>Poland (PL)</td><td>6</td><td>1</td><td>2</td><td>H18, H19</td><td>0.3333 ± 0.2152</td><td>0.0005 ± 0.0007</td><td>- 0.9330</td><td>- 0.0028</td></tr><tr><td>9</td><td>Siberia 2 (SIB2)</td><td>4</td><td>2</td><td>3</td><td>H35—H37</td><td>0.8333 ± 0.2224</td><td>0.0015 ± 0.0015</td><td>- 0.7099</td><td>- 0.8873</td></tr><tr><td>10</td><td>Caucasus Mts. (CAU)</td><td>7</td><td>6</td><td>3</td><td>H20—H22</td><td>0.7143 ± 0.1267</td><td>0.0032 ± 0.0023</td><td>- 0.6351</td><td>1.4625</td></tr><tr><td>11</td><td>North America (NA)</td><td>58</td><td>6</td><td>7</td><td>H1—H7</td><td>0.2281 ± 0.0734</td><td>0.0005 ± 0.0005</td><td>- 1.8875**</td><td>- 6.7157**</td></tr><tr><td colspan="2">Total</td><td>223</td><td>47</td><td>46</td><td>H1—H46</td><td>0.9159 ± 0.0114</td><td>0.0098 ± 0.0052</td><td>- 0.4794</td><td>- 13.864**</td></tr><tr><td colspan="10"><em>Boloria eunomia</em> nDNA (ArgKin)</td></tr><tr><td>1</td><td>W Palaearctic (WP)</td><td>38</td><td>2</td><td>4</td><td>H4—H7</td><td>0.4737 ± 0.0728</td><td>0.0009 ± 0.0009</td><td>0.2595</td><td>- 0.8485</td></tr><tr><td>2</td><td>Scandinavia & Siberia (SC-SIB)</td><td>22</td><td>18</td><td>14</td><td>H8—H21</td><td>0.9177 ± 0.0455</td><td>0.0057 ± 0.0034</td><td>- 1.1663</td><td>- 6.7518**</td></tr><tr><td>3</td><td>North America (NA)</td><td>17</td><td>5</td><td>3</td><td>H1—H3</td><td>0.2279 ± 0.1295</td><td>0.0010 ± 0.0010</td><td>- 1.9433**</td><td>- 0.0292</td></tr><tr><td colspan="2">Total</td><td>77</td><td>30</td><td>21</td><td>H1—H21</td><td>0.8312 ± 0.0299</td><td>0.0085 ± 0.0046</td><td>- 0.5251</td><td>- 3.4341</td></tr><tr><td colspan="10"><em>Boloria selene</em> mtDNA (COI)</td></tr><tr><td>1</td><td>Eurasia (E)</td><td>82</td><td>52</td><td>25</td><td>H9—H19, H21—H34</td><td>0.7064 ± 0.0560</td><td>0.0042 ± 0.0026</td><td>- 2.4878**</td><td>- 16.084**</td></tr><tr><td>2</td><td>Belarus (BR)</td><td>1</td><td>0</td><td>1</td><td>H20</td><td>n.a.</td><td>n.a.</td><td>n.a.</td><td>n.a.</td></tr><tr><td>3</td><td>Spain (SP)</td><td>1</td><td>0</td><td>1</td><td>H18</td><td>n.a.</td><td>n.a.</td><td>n.a.</td><td>n.a.</td></tr><tr><td>4</td><td>North America (NA)</td><td>17</td><td>4</td><td>4</td><td>H1—H4</td><td>0.5000 ± 0.1355</td><td>0.0015 ± 0.0012</td><td>- 0.8161</td><td>- 0.5062</td></tr><tr><td>5</td><td>Alberta (AB)</td><td>2</td><td>1</td><td>2</td><td>H5, H6</td><td>1.0000 ± 0.5000</td><td>0.0017 ± 0.0024</td><td>0.0001</td><td>0.0001</td></tr><tr><td>6</td><td>Wisconsin (W)</td><td>3</td><td>1</td><td>2</td><td>H7, H8</td><td>0.6667 ± 0.3143</td><td>0.0012 ± 0.0014</td><td>0.0001</td><td>0.2007</td></tr><tr><td colspan="2">Total</td><td>106</td><td>84</td><td>34</td><td>H1—H34</td><td>0.8045 ± 0.0376</td><td>0.0145 ± 0.0075</td><td>- 1.4995*</td><td>- 5.9885</td></tr></tbody></table>

Population groups within their Holarctic range are based on SAMOVA analyses.

(N) number of individuals sampled, (S) number of polymorphic sites, (h) number of haplotypes, (Hd) haplotype diversity, (π) nucleotide diversity, with standard deviation (SD), (D) and (Fs) statistic for neutrality test,

** p < 0.01,

* p < 0.05

https://doi.org/10.1371/journal.pone.0214483.t002

Mya, however, relationships among them have to be considered with caution because of the weak branch support.

## *Boloria eunomia* arginine kinase

We obtained 577 bp ArgKin sequences (77 individuals from 47 localities, 21 haplotypes; Fig 3A–3C). The haplotypes from the Nearctic region (H1-3) were again separated from populations of Eurasia. This nuclear marker also showed a strong relation of Scandinavian and Siberian haplotypes (H8-21). However, the general phylogeographic structuring across the Palaearctic region was weaker when compared to COI. In SAMOVA, the highest variance increase of Φ<sub>CT</sub> was for K = 3 (Φ<sub>CT</sub> = 0.65), separating: (1) NA: North America; (2) SC-SIB: Scandinavia + Poland + Siberia and (3) WP: Western Palaearctic (Spain, France, Belgium, Alps, CZ, Poland, Belarus, S Finland, Bulgaria, Caucasus and Ural Mts.) (Fig 2B and Figure B
<figure><img src="image_6.png" /><figcaption>**Fig 3.** Phylogeographic patterns of *Boloria eunomia* based on nuclear (ArgKin) marker. Chronogram calculated in BEAST (A); node bars display 95% highest posterior density interval of the molecular clock and branch labels show posterior probabilities. TCS haplotype network illustrating relationships of the 21 ArgKin haplotypes (B); mutations are shown as hatch marks; the colours correspond to the haplotype distribution map (C).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0214483.g003

in S2 Fig). Overall haplotype diversity was again high (Hd = 0.83) with the highest value for the SC-SIB group (Hd = 0.92). Significant negative values of Tajima's D for NA and Fu's Fs for SC-SIB were observed (Table 2). The BEAST tree topology (Fig 3A) was consistent with both haplotype network (Fig 3B) and SAMOVA (Fig 2B), dividing haplotypes into three groups.

### Boloria selene cytochrome c oxidase

In *B. selene*, we obtained 582 bp COI sequences for 106 individuals from 62 localities and identified 34 haplotypes (Fig 4A and 4B). Haplotype network showed considerable difference between haplotypes from Nearctic and Palaearctic regions (Fig 4B). Haplotypes within the North America were divided into three groups; the first one (H1-4) was characterised by one widespread haplotype occurring throughout Nearctic (H1), the second included samples from Alberta (H5-6) and the last involved Wisconsin haplotypes (H7-8). In contrast, no clear geographic pattern was visible across the Palaearctic region, characterised by star-like network structure (H9-34). In SAMOVA, large increases in $\Phi_{CT}$ occurred at K = 4, 6, 11 and 15 (Figure C in S2 Fig). We evaluated as the best K = 6 ($\Phi_{CT}$ = 0.82) because higher levels of K become less biologically informative (i.e. clusters formed by a single sample). Subsequent population grouping is depicted in Fig 2C: (1) E: Eurasia; (2) BR: Belarus; (3) SP: Spain; (4) NA: North America; (5) W: Wisconsin; (6) AB: Alberta. Total haplotype diversity was again relatively high (Hd = 0.81). In this region, significant negative values in both Tajima's D and Fu's Fs were observed (Table 2).

According to the BEAST chronogram (Fig 4C), COI haplotypes of *B. selene* showed strongly supported division between the Nearctic and Palaearctic populations, which differentiated during the Pleistocene (~2.5 Mya). Populations from the Palaearctic region mostly diverged ~0.5 Mya. The topology of the tree is congruent with the haplotype network (Fig 4B).

### Species distribution modelling (SDM)

**Boloria eunomia.** The best RM was according to AICc = 4.5 and the most parsimonious function of FC was a hinge feature (H) (S2 Table). The average area under the curve (AUC) over the ten replicate runs was 0.870 (± 0.029 SD), indicating a good fit for the model. The selected most important BIOCLIM variables after the jackknife procedure, and the respective
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Phylogeographic patterns of *Boloria selene* based on mitochondrial (COI) marker. Haplotypes distribution on the map (A). TCS haplotype network illustrating relationships of the 34 (COI) haplotypes (B); mutations are shown as hatch marks; the colours correspond to the map. Chronogram calculated in BEAST (C); node bars display 95% highest posterior density interval of the molecular clock and branch labels show posterior probabilities.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0214483.g004
<figure><img src="image_8.png" /><figcaption>**Fig 5.** Species distribution models. SDMs showing the climatically suitable areas for *Boloria eunomia* at (A) present conditions and (B) LGM (about 22,000 years ago) and *Boloria selene* (C) present conditions, (D) LGM. The blue colours show the suitability from 0.14 for *B. eunomia* and 0.064 for *B. selene* (the average minimal training presence threshold, including the least suitable training record) and the red colours from 0.347 for *B. eunomia* and 0.338 for *B. selene* (the average ten percentile training presence threshold, suitability of 90% of the training data). The LGM glacial sheet was redrawn from Ehlers *et al.* [66]; the drawing is not identical to the original shape file and therefore, it is only for illustrative purposes.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0214483.g005

percentages of explained variation, were as follows: BIO5 (maximum temperature of warmest month), 95%; BIO7 (temperature annual range), 3.6%; BIO18 (precipitation of warmest quarter), 1.4%. Based on the selection of variables, the distribution of *B. eunomia* appears to be controlled mostly by temperature. The potential present distribution (Fig 5A) of the species covers
large parts of boreal and Arctic Eurasia and North America, whereas the southern parts of the distribution are modelled as a series of isolated refugia. Some of them are actually occupied by the species (for instance the Alps, Pyrenees, Stara Planina, Caucasus, Altai, Rocky Mts.), whereas some are unoccupied (e.g., Iceland, UK, Sudeten Mts., Carpathians, Pamir Mts.). During the LGM, the model showed widespread distribution across the entire species distribution range reaching further south than in present (Fig 5B).

**Boloria selene.** The selected RM was 4 and the FC was a linear, quadratic, hinge feature (LQH) (S2 Table). The average AUC was 0.793 (± 0.034 SD). The selected BIOCLIM variables, and the respective percentages of explained variation, were: BIO10 (mean temperature of warmest quarter), 50%; BIO18 (precipitation of warmest quarter), 30%; BIO5 (maximum temperature of warmest month), 10.8%; BIO15 (precipitation seasonality), 8.3% and BIO2 (mean diurnal range), 0.9%. The distribution of *B. selene* also appears to be controlled by temperature. The potential present distribution of *B. selene* is similar to *B. eunomia*, but less isolated in southern areas. For instance, occupied isolated distribution is shown for Iberian chain or Caucasus, unoccupied but climatically suitable areas for Iceland or in Central Asian mountains (Fig 5C). The LGM model again showed continuous distribution all across the Holarctic region but with retreat from the northernmost parts of the actual species distribution towards south (Fig 5D).

Discussion

Based on mitochondrial and nuclear markers and species distribution models, we assessed the biogeographic patterns of two Holarctic butterfly species *Boloria eunomia* and *Boloria selene*, providing a comprehensive perspective on the processes affecting these butterflies through the Quaternary.

Palaearctic origin and colonization of America

Central Asia is believed to be the evolutionary cradle for many cold-adapted butterflies such as the genera *Parnassius* Latreille, 1804 [67], *Oeneis* Hübner, [1819] [68], *Erebia* Dalman, 1816 [69], *Coenonympha* Hübner, [1819] [70,71], including *Boloria* Moore, 1900 [49,72]. Our study species, *B. eunomia* and *B. selene* probably also originated in the Palaearctic region [49]. In the case of *B. eunomia*, our results showed the highest genetic diversity within Siberia, supporting the Asian origin (Fig 1B and Table 2). Consequently, both butterflies expanded from Siberia across the other parts of the Holarctic region (Fig 6).

*Boloria selene* most likely colonised North America around 2.5 Mya (Fig 4C) [49] through the second Beringian land bridge where a continuous forest belt connected both continents which might have allowed the dispersal of boreal taxa [73]. This route appears to have been used by other butterflies, with contrasting climatic and habitat affinities, for instance the temperate-arboreal genus *Limenitis* Fabricius, 1807 [74], the warm grasslands generalist *Junonia* Hübner, [1819] [75], the cold-adapted *Oeneis* [66] and other members of the Satyrinae sub-family [69,76,77]. During the further Quaternary cooling and aridification, the third Beringian land bridge re-emerged [73], enabling the colonisation of the Nearctic region by the more cold-tolerant *B. eunomia* around 1.6 Mya (Fig 1C) [49]. Both studied butterflies colonised the Nearctic region only once, long before LGM, they share the significant genetic divergence between populations from Eurasia and North America (Figs 1–4), and were probably isolated during multiple upper Quaternary glacial cycles. The dry climate and tundra habitat might have acted as a barrier to species like *B. eunomia* and *B. selene*, preferring more humid or mesic conditions, respectively. Similar differentiation was observed in several other Lepidoptera (cf. [49,78]) and small mammal species [79–81].
<figure><img src="image_9.png" /><figcaption>**Fig 6.** Hypothetical scenarios of colonization pathways for *Boloria eunomia* (A) and *Boloria selene* (B). Different colours and numbers indicate movement of each haplotype group illustrated in Fig 1 (#1 NA; #2 CS-SIB2; #3 NS-SIB1; #4 A1-CZ, A2, CA-PY, EE-UR, CAU; #5 AR-NI) and Fig 4 (#1 NA; #2 W; #3 E).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0214483.g006
## Biogeography in North America

Within North America, both species differ in their genetic structures. Our results imply only one glacial refugium for *B. eunomia*. Star-like pattern in haplotype network (Fig 1B), low genetic diversity across vast area of the Nearctic and non-significant neutrality tests (Table 2) suggest recent population expansion, probably from never glaciated Alaska [82] or southern Canada and northern USA, south of the Laurentide Ice Sheet. Unfortunately, without samples from southern parts of Rocky Mts., we are unable to choose between these alternatives.

In comparison with *B. eunomia*, *B. selene* colonised North America much earlier and is characterised by three different lineages within the region (Figs 2C and 4B), implying glacial survival not only in Alaska but also south of the Laurentide ice-sheet in the so-called 'Driftless Region' and in Rocky Mts. Previous phylogeographic and palaeontological studies suggest that the 'Driftless Region' (south-western Wisconsin, south-eastern Minnesota and north-eastern Iowa) served as a refugium for many other plant and animal species [83–90]. This is supported by the occurrence of unique haplotypes in the Wisconsin *B. selene* population (Fig 4B) and the highest genetic diversity within the species' Nearctic range (Table 2). The third *B. selene* lineage is formed by two private haplotypes from Alberta (Fig 4B) suggesting another refugium and/or physical barriers within Rocky Mts. The split between haplotypes from Alberta and British Columbia was also observed in two moth species *Malacosoma disstria* (Hübner, 1820) and *M. californica* (Packard, 1864) [91]. All these possible LGM refugia are in concordance with the SDM models (Fig 5B and 5D) showing suitable areas within Beringia, north-west Rocky Mts. and south of the Laurentide ice-shield.

## Boloria eunomia—Palaearctic phylogeography

Contrary to the findings of Yakovlev et al. [92], our results confirm a high degree of phylogeographic structure within the *Boloria eunomia* Eurasian range. A high within-range genetic differentiation can be observed in other habitat specialist butterflies, for instance the semi-migratory species of grasslands *Hesperia comma* (Linnaeus, 1758) [93], the open-woodland specialist *Parnassius mnemosyne* (Linnaeus, 1758) [94], the polycentric dry grassland species *Melitaea cinxia* (Linnaeus, 1758) [95] or the steppic *Proterebia afra* (Fabricius, 1787) [96]. We identified ten population groups of *B. eunomia* (Fig 2A) within the Palaearctic region, mostly differentiating around 0.5 Mya (Figs 1C and 3A), corresponding with the onset of the long 100 kya glacial cycles of the Middle and Late Pleistocene. Only two lineages spread over the vast area of Siberia, leaving the highest geographic diversity to Europe, indicating different species' histories across continents.

**Dynamics in Siberia, Fennoscandia and Eastern Europe.** Two extensive lineages exist in Siberia (NS-SIB1, SIB2; Figs 1A–1C and 2A), harbouring the highest haplotype diversity within the species (Table 2), possibly as a result of an interconnected network of populations during both glacials and interglacials. Siberia endured much less severe glaciations compared to Europe [97]. Moreover, this scenario is supported by the wide potential distribution during the LGM (Fig 5B), and the lack of geography-dependent genetic differentiation. Also, the archaeobotanical data and pollen records confirm high abundances of *B. eunomia* larval food plants (*Polygonum bistorta* and *P. viviparum*) in the European tundra about 15 kya (e.g. [98]). Such pattern had been already reported in another boreo-montane butterfly *Lycaena helle* [13].

As shown in other boreo-montane species [14,99–101], the two eastern lineages are responsible for post-glacial (re)colonisation of Fennoscandia (Figs 1 and 6A). The third lineage (EE-UR) reaching Scandinavia encompasses large parts of Eastern Europe with affinities to Ural Mts. and the Balkans (Stara Planina Mts.) (Figs 1A–1C and 2A). The proximity of
Bulgarian and Polish populations may suggest a hypothetical existence of now extinct Carpa-
thian populations of this species.

Not only Scandinavia but also Poland might have acted as contact zone for various lineages of *B. eunomia* as hypothesised by Krzywicki [102]. Besides of the EE-UR lineage, it also consists of an independent Polish clade (PL) and it shares haplotypes with NS-SIB1 lineage (Figs 1 and 2A and 2B and Fig B in S1 Fig). Similar phenomenon was observed for other butterflies [94,103,104]. It is also worth noting that in Poland two morphological forms of the butterfly showing distinct ecological characteristics co-occurs [105]. Detailed studies on relationships between those populations would be important for understanding of evolution of host races and validity of subspecies division suggested e.g. by Tolman & Lewington [32].

**Diversification among European mountain ranges.** We detected strong differentiation among most European mountain ranges (Cantabrian, Pyrenees, Ardennes, Alps, Bohemian Forest, Stara Planina and Caucasus Mts.) indicating restricted gene flow and long-term isolation of European populations, sufficient to fix population-specific haplotypes (Fig 1B). This genetic structuring contradicts our SDMs showing species' widespread distribution during the LGM. However, more severe glaciations in the past (Elsterian/Mindel and Saalian/Riss) could have forced even cold-adapted organisms as *B. eunomia* to retreat within several refugia [106]. Mountains offer a heterogeneous topography with favourable microhabitats allowing the permanent persistence of various species over several glacial cycles by altitudinal shifts. At the same time, mountain ranges acted as barrier for the gene flow among populations, resulting in high geographic differentiation. This is the common pattern in the cold-adapted species in general [107] including boreo-montane butterflies [13,14].

Two *B. eunomia* lineages evolved within the Alps (A1-CZ, A2; Figs 1B and 2A), as a result of separation during the range contractions. The eastern lineage A1-CZ colonized the Bohemian Forest (Czech Republic), where the species persisted continuously at least since the end of the last glacial period [25]. Within south-western Europe, the SAMOVA analysis (Fig 2A) recognised two more closely related lineages (Fig 1B and 1C), one including Belgian and French populations (AR-NI), the other covering the Iberian peninsula (CA-PY), confirming the morphometric differentiation of these populations and the lower genetic variation of Ardennes populations compared to Pyrenees' ones [108]. Furthermore, the haplotype network and Bayesian tree (Fig 1B and 1C) detected a differentiation between Cantabrian Mts. and Pyrenees, another example of population separation among different mountain ranges, reported for many plant and animal species [11,109–112]. Population in Caucasus Mts. (CAU) seems to be of the European origin as both mtDNA (Fig 1C) and nDNA (Fig 3D and 3E) suggest affinity to the western rather than the eastern samples. High haplotype diversity (Fig 1B and Table 2) and SDMs (Fig 5A and 5B), showing relative isolation of Caucasus during both present and LGM periods, also supports long-lasting persistence and isolation of the Caucasian population.

***Boloria selene***—**Evolutionary history within Eurasia**

Within the Palaearctic region, both species are characterised by opposite genetic patterns to the ones in the Nearctic. Populations of *Boloria selene* are homogenously distributed and interconnected throughout the large parts of its Eurasian distribution range, reflected by one dominant haplotype (Fig 4B). Shallow mitogenetic/geographic structure in Lepidoptera across large species ranges similar to *B. selene* has been reported for several other butterfly species: two Holarctic—*Parnassius phoebus* (Fabricius, 1793) [113] and *Parasemia plantaginis* (Linnaeus, 1758) [114]; and three Palaearctic species—*Polygonia c-album* (Linnaeus, 1758) [115], *Lopinga achine* (Scopoli, 1763) [116] or *Aglais urticae* (Linnaeus, 1758) [117]. High haplotype diversity, significant negative neutrality tests (Table 2), the wide geographic distribution of common
haplotypes, presence of many private haplotypes, and at the same time the lack of geographic structure (Fig 4A and 4B) suggest ongoing gene flow among *B. selene* populations with a recent population expansion (Fig 6B). This is in agreement with mark-recapture data, as the species easily disperse and demonstrates a low site fidelity [118]. On the other hand, the dispersal of *B. eunomia* is restricted, especially in a fragmented landscape [119]. SDMs (Fig 5C and 5D) indicated widespread distribution in both present and LGM, showing the interconnection of populations through the Quaternary.

## Conclusions

Despite the considerable overlap in their distribution, the genetic patterns in the two *Boloria* species notably differed. *B. eunomia* is characterised by a high genetic structuring across Europe, but panmixia in Asia and a lack of genetic differentiation across the Nearctic region. In contrast, *B. selene* is genetically strongly differentiated across North America and shows gene flow and thus low geographic differentiation across Eurasia. Generalist taxa such as the butterfly *B. selene* tend to have high genetic diversities and a low genetic differentiation among local populations. On the other hand, habitat specialists such as the butterfly species *B. eunomia* are affected from population isolation, which causes strong genetic differentiation [120]. This coherence becomes approved by the genetic signature obtained for the Palaearctic range for *B. selene* and for the European range for *B. eunomia*. Therefore, global evolutionary history of boreo-montane and boreo-temperate butterflies cannot be explained by the single biogeographic scenario proposed for cold-adapted species by Schmitt [8] but rather by combination of them, emphasising the different species' histories in distinct parts of the distribution range. We highlight the importance of studying the genetic diversity within the entire species distribution ranges, which helps us to understand how the evolutionary history of organisms influenced their present distribution, and, thus to prioritize areas of conservation.

## Supporting information

**S1 Fig. Detailed maps of haplotype diversity of *Boloria eunomia* within Europe based on COI.** Two lineages occurring in the Alps and the Czech Republic (A). Haplotype diversity in north-eastern Europe showing the contact zones in Scandinavia and Poland (B).
(TIF)

**S2 Fig. Graphs from SAMOVA analyses.** *Boloria eunomia* cytochrome c oxidase subunit I (COI; 11 populations) (A); and arginine kinase gene (ArgKin; 3 populations) (B); *Boloria selene* cytochrome c oxidase subunit I (COI; 6 populations) (C). K, number of clusters; Φ<sub>CT</sub>, fixation index.
(TIF)

**S1 Table. Overview of samples of both *Boloria eunomia* and *Boloria selene*.** Specimens localities; coordinates; COI, ArgKin and Wingless haplotypes, SAMOVA groups and GenBank accession numbers.
(XLSX)

**S2 Table. MaxEnt.** Comparisons of (A) Regularization Multiplier (RM) parameter and (B) different feature classes (FC) of MaxEnt models by ENMeval R package using AICc, used for calculations of Spatial Distribution Models of *Boloria eunomia* and *Boloria selene*. The lowest AICc values indicate the best models. (C) Relative gain of the MaxEnt model after removing each of the variables. Only variables with positive gain larger than 0.01 were used.
(XLSX)
**Acknowledgments**

We thank the following persons for providing samples and ideas: Bernard Barascud, Joe Belicek, Jiří Beneš, Stoyan Beshkov, Martin Česánek, Oldřich Čížek, Jiří Darebník, Henri Desclimon, Stefan Ehl, Dmitriy Goshko, Vladimír Hula, Dmitriy Morgun, Martin Konvička, Dan Leština, Paul A. Opler, Krzysztof Pałka, Alois Pavlíčko, Anežka Pavlíková, Jure Rekelj, Nikolay Rubin, Juriy Shevnin, Thomas Schmitt, Jiří Skala, Lukáš Spitzer, Pavla Tájková, Derek Sikes, Kath M. Daly. We thank Martin Husemann (Hamburg, Germany) for constructive comments on our manuscript.

**Author Contributions**

**Conceptualization:** Jana Maresova, Zdenek Faltynek Fric.

**Data curation:** Jana Maresova, Jan Christian Habel, Gabriel Neve, Marcin Sielezniew, Agata Kostro-Ambroziak, Zdenek Faltynek Fric.

**Formal analysis:** Jana Maresova, Zdenek Faltynek Fric.

**Funding acquisition:** Zdenek Faltynek Fric.

**Investigation:** Jana Maresova, Jan Christian Habel, Gabriel Neve, Marcin Sielezniew, Alena Bartonova, Agata Kostro-Ambroziak, Zdenek Faltynek Fric.

**Methodology:** Jana Maresova, Jan Christian Habel, Gabriel Neve, Marcin Sielezniew, Zdenek Faltynek Fric.

**Project administration:** Jana Maresova.

**Resources:** Jana Maresova, Jan Christian Habel, Gabriel Neve, Marcin Sielezniew, Alena Bartonova, Agata Kostro-Ambroziak, Zdenek Faltynek Fric.

**Software:** Jana Maresova, Zdenek Faltynek Fric.

**Supervision:** Zdenek Faltynek Fric.

**Validation:** Jana Maresova, Jan Christian Habel, Gabriel Neve, Marcin Sielezniew, Alena Bartonova, Zdenek Faltynek Fric.

**Visualization:** Jana Maresova.

**Writing – original draft:** Jana Maresova.

**Writing – review & editing:** Jana Maresova, Jan Christian Habel, Gabriel Neve, Marcin Sielezniew, Alena Bartonova, Zdenek Faltynek Fric.

**References**

1. Avise JC. Phylogeography—the history and formation of species. Cambridge, Massachusetts, London: Harvard University Press; 2000.

2. Hewitt GM. The structure of biodiversity—insights from molecular phylogeography. Front Zool. 2004; 1:4. https://doi.org/10.1186/1742-9994-1-4 PMID: 15679920

3. Hewitt GM. Post-glacial re-colonization of European biota. Biol J Linn Soc Lond. 1999; 68:87–112.

4. Seddon JM, Santucci F, Reeve NJ, Hewitt GM. DNA footprints of European hedgehogs, *Erinaceus europaeus* and *E. concolor*: Pleistocene refugia, postglacial expansion and colonization routes. Mol Ecol. 2001; 10:2187–2198. PMID: 11555261

5. Varga Z, Schmitt T. Types of areal and oretundral disjunctions in the western Palearctic. Biol J Linn Soc Lond. 2008; 93:415–430.

6. Ronikier M, Schneeweiss GM, Schoönswetter P. The extreme disjunction between Beringia and Europe in *Ranunculus glacialis* s. l. (Ranunculaceae) does not coincide with the deepest genetic split
—a story of the importance of temperate mountain ranges in arctic-alpine phylogeography. Mol Ecol. 2012; 21:5561–5578. https://doi.org/10.1111/mec.12030 PMID: 23057869

7. Schmitt T, Varga Z. Extra-Mediterranean refugia: The rule and not the exception? Front Zool. 2012; 9:22. https://doi.org/10.1186/1742-9994-9-22 PMID: 22953783

8. de Lattin G. Grundriss der Zoogeographie. Fischer, Stuttgart: Fischer; 1967.

9. Schmitt T. Molecular biogeography of Europe: Pleistocene cycles and postglacial trends. Front Zool. 2007; 4:11. https://doi.org/10.1186/1742-9994-4-11 PMID: 17439649

10. Després L, Loriot S, Gaudeul M. Geographic pattern of genetic variation in the European globeflower *Trollius europaeus* L. (Ranunculaceae) inferred from amplified fragment length polymorphism markers. Mol Ecol. 2002; 11:2337–2347. PMID: 12406244

11. Kramp K, Huck S, Niketić M, Tomović G, Schmitt T. Multiple glacial refugia and complex postglacial range shifts of the obligatory woodland plant *Polygonatum verticillatum* (Convallariaceae). Plant Biol (Stuttg). 2009; 11:392–404.

12. Drees C, Husemann M, Homburg K, Brandt P, Dieker P, Habel JC, et al. Molecular analyses and species distribution models indicate cryptic northern mountain refugia for a forest-dwelling ground beetle. J Biogeogr. 2016; 43:2223–2236.

13. Habel JCH, Schmitt T, Meyer M, Finger A, Rödder D, Assmann T, et al. Biogeography meets conservation: the genetic structure of the endangered lycaenid butterfly *Lycaena helle* (Denis & Schiffermüller, 1775). Biol J Linn Soc Lond. 2010; 101:155–168.

14. Kramp K, Cizek O, Madeira PM, Ramos AA, Konvicka M, Castilho R, et al. Genetic implications of phylogeographical patterns in the conservation of the boreal wetland butterfly *Colias palaeno* (Pieridae). Biol J Linn Soc Lond. 2016; 119:1068–1081.

15. Nève G, Barascud B, Hughes R, Aubert J, Descimon H, Lebrun P, et al. Dispersal, colonization power and metapopulation structure in the vulnerable butterfly *Proclossiana eunomia* (Lepidoptera, Nymphalidae). J Appl Ecol. 1996; 33:14–22.

16. Binney HA, Willis KJ, Edwards ME, Bhagwat SA, Anderson PM, Andreev AA, et al. The distribution of late-Quaternary woody taxa in northern Eurasia: evidence from a new macrofossil database. Quat Sci Rev. 2009; 28:2445–2464.

17. Ehlers J, Gibbard PL. The extent and chronology of Cenozoic global glaciation. Quat Int 2007; 164:6–20.

18. Tarasov PE, Volkova VS, Webb T III, Guiot J, Andreev AA, Bezusko LG, et al. Last glacial maximum biomes reconstructed from pollen and plant macrofossil data from northern Eurasia. J Biogeogr. 2000; 27:609–620.

19. Abbott RJ, Brochmann C. History and evolution of the Arctic flora: in the footsteps of Eric Hultén. Mol Ecol. 2003; 12:299–313. PMID: 12535083

20. Goetcheus VG, Birks HH. Full-glacial upland tundra vegetation preserved under tephra in the Berin-gian National Park, Seward Peninsula, Alaska. Quat Sci Rev. 2001; 20:135–147.

21. Alsos IG, Engelskjøn T, Gielly L, Taberlet P, Brochmann C. Impact of ice ages on circumpolar molecular diversity: insights from an ecological key species. Mol Ecol. 2005; 14:2739–2753. https://doi.org/10.1111/j.1365-294X.2005.02621.x PMID: 16029475

22. Eidesen PB, Alsos IG, Popp M, Stensrud O, Suda J, Brochmann C. Nuclear vs. plastid data: complex Pleistocene history of a circumpolar key species. Mol Ecol. 2007; 16:3902–3925. https://doi.org/10.1111/j.1365-294X.2007.03425.x PMID: 17850553

23. Mardulyn P, Mikhailov Y, Pasteels JM. Testing phylogeographic hypotheses in a Euro-Siberian cold-adapted leaf beetle with coalescent simulations. Evolution. 2009; 63:2717–2729. https://doi.org/10.1111/j.1558-5646.2009.00755.x PMID: 19545269

24. Nève G, Barascud B, Descimon H, Baguette M. Gene flow rise with habitat fragmentation in the bog fritillary butterfly (Lepidoptera: Nymphalidae). BMC Evol Biol. 2008; 8:4. https://doi.org/10.1186/1471-2148-8-84 PMID: 18366652

25. Nève G, Pavlíčko A, Konvička M. Loss of genetic diversity through spontaneous colonization in the bog fritillary butterfly, *Proclossiana eunomia* (Lepidoptera: Nymphalidae) in the Czech Republic. Eur J Entomol. 2009; 106:11–19.

26. Tuzov VK, Bozano GC. Guide to the butterflies of the Palearctic region. Nymphalidae, Part II, Tribe Argynnini: *Boloria, Proclossiana, Glossiana*. Milano: Omnes Artes; 2006.

27. Petersen B. Some trends of speciation in the cold-adapted Holarctic fauna. Zoologiska Bidrag från Uppsala. 1956; 30:233–314.

28. Jakšić P, Van Swaay C, urić M. *Boloria eunomia* (Esper, 1799): a new species for Serbia (Nymphalidae). Nota Lepidopterol. 2007; 30:64–70.
29. Kudrna O, Pennerstorfer J, Lux K. Distribution atlas of European butterflies and skippers. 3rd ed. Schwanfeld: Wissenschaftlicher Verlag Peks; 2015.

30. Schtickzelle N, Turlure C, Baguette M. Grazing management impacts on the viability of the threatened bog fritillary butterfly *Proclossiana eunomia*. *Biol Conserv.* 2007; 136:651–660.

31. Klimczuk P, Sielezniew M. Unexpected differences in butterfly diversity between two peat bogs in the same area. *Pol J Entomol.* 2017; 108:251–273.

32. Tolman T, Lewington R. Collins butterfly guide: the most complete field guide to the butterflies of Britain and Europe. 3rd ed. London: HarperCollins UK; 2008.

33. Chernov Yul, Tatarinov AG. Butterflies (Lepidoptera, Rhopalocera) in the Arctic fauna. *Entomol Rev.* 2006; 86(7):760–786.

34. Beneš J, Konvika M, Dvořák J, Fric Z, Havelda Z, Pavlíčko A, et al. Butterflies of the Czech Republic: Distribution and Conservation I, II. Prague: SOM; 2002.

35. Layberry RA, Hall PW, Lafontaine J. The Butterflies of Canada. Toronto: University of Toronto Press; 1998.

36. Tshikolovets VV. Butterflies of Europe & the Mediterranean area. Kiev: Tshikolovets; 2011.

37. Wahlberg N, Wheat CW. Genomic outposts serve the phylogenomic pioneers: designing novel nuclear markers for genomic DNA extractions of Lepidoptera. *Syst Biol.* 2008; 57:231–242. https://doi.org/10.1080/10635150802033006 PMID: 18398768

38. Monteiro A, Pierce NE. Phylogeny of *Bicyclus* (Lepidoptera: Nymphalidae) inferred from COI, COII, and EF-1alpha gene sequences. *Mol Phylogenet Evol.* 2001; 18:264–281. https://doi.org/10.1006/mpev.2000.0872 PMID: 11161761

39. Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. *Bioinformatics.* 2012; 28:1647–1649. https://doi.org/10.1093/bioinformatics/bts199 PMID: 22543367

40. Clement M, Posada D, Crandall K. TCS: a computer program to estimate gene genealogies. *Mol Ecol.* 2000; 9:1657–1660. PMID: 11050560

41. Leigh JW, Bryant D. PopART: full-feature software for haplotype network construction. *Methods Ecol Evol.* 2015; 6:1110–1116.

42. Dellicour S, Mardulyn P. SPADS 1.0: a toolbox to perform spatial analyses on DNA sequence data sets. *Mol Ecol Resour.* 2014; 14:647–651. https://doi.org/10.1111/1755-0998.12200 PMID: 24215429

43. Excoffier L, Laval G, Schneider S. Arlequin vers. 3.0: an integrated software package for population genetics data analysis. *Evol Bioinform Online.* 2005; 1:47–50.

44. Excoffier L, Lischer HEL. Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. *Mol Ecol Resour.* 2010; 10:564–567. https://doi.org/10.1111/j.1755-0998.2010.02847.x PMID: 21565059

45. Lanfear R, Calcott B, Ho Sy, Guindon S. Partition-finder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. *Mol Biol Evol*, 2012, 29:1695–1701. https://doi.org/10.1093/molbev/mss020 PMID: 22319168

46. Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. *Mol Biol Evol.* 2012; 29:1969–1973. https://doi.org/10.1093/molbev/mss075 PMID: 22367748

47. Wahlberg N. That awkward age for butterflies: insights from the age of the butterfly subfamily Nymphalinae. *Syst Biol.* 2006; 55:703–714. https://doi.org/10.1080/10635150600913235 PMID: 16952908

48. Wahlberg N, Wheat CW, Peña C. Timing and patterns in the taxonomic diversification of Lepidoptera (butterflies and moths). *PLoS One.* 2013; 8(11):e80875.

49. Simonsen TJ, Wahlberg N, Warren AD, Sperling FAH. The evolutionary history of *Boloria* (Lepidoptera: Nymphalidae): phylogeny, zoogeography and larval foodplant relationships. *Syst Biodivers.* 2010; 8:513–529.

50. Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. *BMC Evol Biol.* 2007; 7:214. https://doi.org/10.1186/1471-2148-7-214 PMID: 17996036

51. Miller MA, Pfeiffer W, Schwartz T. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. In: *Proceedings of the Gateway Computing Environments Workshop (GCE). New Orleans: IEEE; 2010.*

52. Tshikolovets VV, Bidzilya VO, Golovushkin MI. The butterflies of Transbaikal Siberia. Kiev, Pardubice: Tshikolovets; 2002.

53. Tshikolovets VV, Yakovlev RV, Bálint Z. The butterflies of Mongolia. Kiev, Pardubice: Tshikolovets; 2009.
54. Tshikolovets VV, Yakovlev RV, Kosterin OE. The butterflies of Altai, Sayans and Tuva (S.-W. Siberia). Kiev, Pardubice: Tshikolovets; 2009.

55. Tshikolovets VV, Nekrutenko YP. The butterflies of Caucasus and Transcaucasia (Armenia, Azerbaijan, Georgia and Russia). Kiev, Pardubice: Tshikolovets; 2012.

56. Tshikolovets VV, Kosterin O, Gorbunov P, Yakovlev R. The Butterflies of Kazakhstan. Kiev, Pardubice: Tshikolovets; 2016.

57. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005; 25:1965–1978.

58. Phillips SJ, Dudík M, Schapire RE. [Internet] Maxent software for modeling species niches and distributions (Version 3.4.1). [cited 2018 Nov 29]. http://biodiversityinformatics.amnh.org/open_source/maxent/.

59. Hijmans RJ, Phillips S, Leathwick J, Elith J. dismo: Species Distribution Modeling. R package version 1.1–4. 2017. https://cran.r-project.org/web/packages/dismo/dismo.pdf

60. Hijmans RJ. raster: Geographic Data Analysis and Modeling. R package version 2.8–4. 2018. https://cran.r-project.org/web/packages/raster/raster.pdf

61. Fourcade Y, Engler JO, Rodder D, Secondi J. Mapping species distributions with MAXENT using a geographically biased sample of presence data: A performance assessment of methods for correcting sampling bias. PLoS One. 2014; 9:e97122. https://doi.org/10.1371/journal.pone.0097122 PMID: 24818607

62. QGIS Development Team. [Internet] QGIS Geographic information system. Open source geospatial foundation project. [cited 2018 Nov 29]. http://qgis.osgeo.org.

63. Muscarella R, Galante PJ, Soley-Guardia M, Boria RA, Kass JM, Uriarte M, et al. ENMeval: an R package for conducting spatially independent evaluations and estimating optimal model complexity for maxent eco-logical niche models. Methods Ecol Evol. 2014; 5:1198–1205.

64. Fu Y-X. Statistical tests of neutrality of mutations against population growth, hitch-hiking, and background selection. Genetics. 1997; 147:915–925. PMID: 9335623

65. Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989; 123:585–595. PMID: 2513255

66. Ehlers J, Gibbard PL, Hughes PD. Quaternary glaciations—extent and chronology. Amsterdam: Elsevier; 2011.

67. Todisco V, Gratton P, Cesaroni D, Sbordoni V. Phylogeography of *Parnassius apollo*: Hhints on taxonomy and conservation of a vulnerable glacial butterfly invader. Biol J Linn Soc Lond. 2010; 101:169–183.

68. Kleckova I, Cesanek M, Fric Z, Pellissier L. Diversification of the cold-adapted butterfly genus *Oeneis* related to Holarctic biogeography and climatic niche shifts. Mol Phylogenet Evol. 2015; 92:255–265. https://doi.org/10.1016/j.ympev.2015.06.012 PMID: 26166775

69. Peña C, Witthauer H, Klečková I, Fric Z, Wahlberg N. Adaptive radiations in butterflies: evolutionary history of the genus *Erebia* (Nymphalidae: Satyrinae). Biol J Linn Soc Lond. 2015; 116:449–467.

70. Kodandaramaiah U, Wahlberg N. Phylogeny and biogeography of *Coenonympha* butterflies (Nymphalidae: Satyrinae)—patterns of colonization in the Holarctic. Syst Entomol. 2009; 34:315–323.

71. Kodandaramaiah U, Braby MF, Grund R, Müller CJ, Wahlberg N. Phylogenetic relationships, biogeography and diversification of Coenonymphina butterflies (Nymphalidae: Satyrinae): intercontinental dispersal of a southern Gondwanan group?. Syst Entomol. 2018; 43:798–809.

72. Varga Z. Das Prinzip der areal-analytischen Methode in der Zoogeographie und die Faunenelemente-Einteilung der europäischen Tagschmetterlinge (Lepidoptera: Diurna). Acta Biologica Debrecina. 1977; 14:223–285.

73. Sanmartín I, Enghoff H, Ronquist F. Patterns of animal dispersal, vicariance and diversification in the Holarctic. Biol J Linn Soc Lond. 2001; 73:345–390.

74. Mullen SP. Wing pattern evolution and the origins of mimicry among North American admiral butterflies (Nymphalidae: Limenitis). Mol Phylogenet Evol. 2006; 39:747–758. https://doi.org/10.1016/j.ympev.2006.01.021 PMID: 16500119

75. Kodandaramaiah U, Wahlberg N. Out-of-Africa origin and dispersal-mediated diversification of the butterfly genus *Junonia* (Nymphalidae: Nymphalinae). J Evol Biol. 2007; 20:2181–2191. https://doi.org/10.1111/j.1420-9101.2007.01425.x PMID: 17887973

76. Peña C, Nylin S, Freitas AVL, Wahlberg N. Biogeographic history of the butterfly subtribe Euptychiina (Lepidoptera, Nymphalidae, Satyrinae). Zool Scr. 2010; 39:243–258.

77. Peña C, Nylin S, Wahlberg N. The radiation of Satyrini butterflies (Nymphalidae: Satyrinae): a challenge for phylogenetic methods. Zool J Linn Soc. 2011; 161:64–87.
78. Lafontaine JD, Wood DM. A zoogeographic analysis of the Noctuidae (Lepidoptera) of Beringia, and some inferences about past Beringian habitat. *Memoirs of the Entomological Society of Canada*. 1988; 144:109-123.

79. Galbreath KE, Cook JA. Genetic consequences of Pleistocene glaciations for the tundra vole (Microtus oeconomus) in Beringia. Mol Ecol. 2004; 13:135-148. PMID: 14653795

80. Hope AG, Takebayashi N, Galbreath KE, Talbot SL, Cook JA, Riddle B. Temporal, spatial and ecological dynamics of speciation among amphi-Beringian small mammals. J Biogeogr. 2013; 40:415-429.

81. Kohli BA, Fedorov VB, Waltari E, Cook JA. Phylogeography of a Holarctic rodent (*Myodes rutilus*): Testing high-latitude biogeographical hypotheses and the dynamics of range shifts. J Biogeogr. 2015; 42:377-389.

82. Dyke AS. (2004) An outline of North American Deglaciation with emphasis on central and northern Canada. In: Ehlers J, Gibbard PL, editors. Quaternary Glaciations—Extent and Chronology, Part II. Amsterdam: Elsevier; 2004. p. 373-424.

83. Beatty GE, Provan JIM. Refugial persistence and postglacial recolonization of North America by the cold-tolerant herbaceous plant *Orthilia secunda*. Mol Ecol. 2010; 19:5009-21. https://doi.org/10.1111/j.1365-294X.2010.04859.x PMID: 20958819

84. Godbout J, Jaramillo-Correa JP, Beaulieu J, Bosquet J. A mitochondrial DNA minisatellite reveals the postglacial history of jack pine (*Pinus banksiana*), a broad-range North American conifer. Mol Ecol. 2005; 14:3497-3512. https://doi.org/10.1111/j.1365-294X.2005.02674.x PMID: 16156818

85. Jackson ST, Webb RS, Anderson KH, Overpeck JT, Webb T III, Williams JW, et al. Vegetation and environment in Eastern North America during the Last Glacial Maximum. Quat Sci Rev. 2000; 19:489-508.

86. Jaramillo-Correa JP, Beaulieu J, Bousquet J. Variation in mitochondrial DNA reveals multiple distant glacial refugia in black spruce (*Picea mariana*), a transcontinental North American conifer. Mol Ecol. 2004; 13:2735-47. https://doi.org/10.1111/j.1365-294X.2004.02258.x PMID: 15315685

87. Latutrie M, Bergeron Y, Tremblay F. Fine-scale assessment of genetic diversity of trembling aspen in northwestern North America. BMC Evol Biol. 2016; 16:231. https://doi.org/10.1186/s12862-016-0810-1 PMID: 27782797

88. Lee-Yaw JA, Irwin JT, Green DM. Postglacial range expansion from northern refugia by the wood frog, *Rana sylvatica*. Mol Ecol. 2008; 17:867-884. https://doi.org/10.1111/j.1365-294X.2007.03611.x PMID: 18179428

89. Li P, Li M, Shi Y, Zhao Y, Wan Y, Fu C, et al. Phylogeography of North American herbaceous *Smilax* (Smilacaceae): combined AFLP and cpDNA data support a northern refugium in the Driftless Area. Am J Bot. 2013; 100:801-814. https://doi.org/10.3732/ajb.1200250 PMID: 23538874

90. Rowe KC, Heske EJ, Brown PW, Paige KN. Surviving the ice: Northern refugia and postglacial recolonization. Proc Natl Acad Sci U S A. 2004; 101:10355-10359. https://doi.org/10.1073/pnas.0401338101 PMID: 15249673

91. Lait LA, Hebert P. Phylogeographic structure in three North American tent caterpillar species (Lepidoptera: Lasiocampidae): *Malacosoma americana*, *M. californica*, and *M. disstria*. PeerJ. 2018; 6: e4479. https://doi.org/10.7717/peerj.4479 PMID: 29576956

92. Yakovlev RV, Shapoval NA, Kuftina GN, Kulak AV, Kovalev SV. Notes on the molecular taxonomy of the *Proclossiana eunomia* complex (Lepidoptera, Nymphalidae: Argynnini): analysis of DNA barcodes. Ukr J Ecol. 2018; 8(1):222-232.

93. Forister ML, Fordyce JA, Shapiro AM. Geological barriers and restricted gene flow in the holarctic skipper *Hesperia comma* (Hesperiidae). Mol Ecol. 2004; 13:3489-3499. https://doi.org/10.1111/j.1365-294X.2004.02349.x PMID: 15488006

94. Gratton P, Konopinski MK, Sbordoni V. Pleistocene evolutionary history of the Clouded Apollo (*Parnassius mnemosyne*): genetic signatures of climate cycles and a 'time-dependent' mitochondrial substitution rate. Mol Ecol. 2008; 17:4248-4262. https://doi.org/10.1111/j.1365-294X.2008.03901.x PMID: 18986502

95. Wahlberg N, Saccheri I. The effects of Pleistocene glaciations on the phylogeography of *Melitaea cinxia* (Lepidoptera: Nymphalidae). Eur J Entomol. 2007; 104:675-684.

96. Bartonova A, Konvicka M, Korb S, Kramp K, Schmitt T, Faltynek Fric Z. (2018) Range dynamics of Palaearctic steppe species under glacial cycles: the phylogeography of *Proterebia afra* (Lepidoptera: Nymphalidae: Satyrinae). Biol J Linn Soc Lond. 2018; 125(4):867-884.

97. Svendsen JI, Alexanderson H, Astakhov VI, Demidov I, Dowdeswell JA, Funder S, et al. Late Quaternary ice sheet history of Eurasia. Quat Sci Rev. 2004; 23:1229-1271.

98. Mägdefrau K. Paläobiologie der Pflanzen. Jena: Verlag Gustav Fischer; 1953.
99. Brunhoff C, Yoccoz NG, Ims RA, Jaarola M. Glacial survival or late glacial colonization? Phylogeography of the root vole (*Microtus oeconomus*) in north-west Norway. *J Biogeogr.* 2006; 33:2136–2144.

100. Mikkola K, Lafontaine JD, Kononenko VS. Zoogeography of the Holarctic species of the Noctuidae (Lepidoptera): importance of the Beringian refuge. *Entomol Fenn.* 1991; 2:157–173.

101. Quinzin MC, Normand S, Dellicour S, Svenning J-C, Mardulyn P. Glacial survival of trophically linked boreal species in northern Europe. *Proc R Soc Lond B Biol Sci.* 2017; 284:20162799.

102. Krzywicki M. Fauna Papilionoidea i Hesperioidea (Lepidoptera) Puszczy Białowieskiej. *Annales Zoologici, Polska Akademia Nauk*. 1967; 25:1–213.

103. Sielezniew M, Ponikwicka-Tyszko D, Ratkiewicz M, Dziekanska I, Kostro-Amnroziak A, Rutkowski R. Divergent patterns in the mitochondrial and nuclear diversity of the specialized butterfly *Plebejus argus* (Lepidoptera: Lycaenidae). *Eur J Entomol.* 2011; 108:537–545.

104. Ratkiewicz M, Jaroszewicz B. Allopatric origins of sympatric forms: the skippers *Carterocephalus palaemon palaemon*, *C. p. tolli* and *C. silvicolus*. *Ann Zool Fennici*. 2006; 43:285–94.

105. Sielezniew M, Kostro-Ambroziak A, Klimczuk P, Deoniziak K, Pałka K, Nowicki P. Habitat-related differences in the adult longevity of two ecotypes of a specialized butterfly. *J Zool*. 2018;1–12.

106. Ehlers J, Astakhov V, Gibbard PL, Mangerud J, Svendsen JI. Middle Pleistocene in Eurasia. In: Elias S, Mock C, editors. *Encyclopedia of Quaternary Science*. 2nd ed. Amsterdam: Elsevier; 2013. p. 172–179.

107. Schmitt T. Biogeographical and evolutionary importance of the European high mountain systems. *Front Zool.* 2009; 6:9. https://doi.org/10.1186/1742-9994-6-9 PMID: 19480666

108. Nève G, Barascud B, Windig JJ. Population biology of *Proclossiana eunomia* (Nymphalidae): Preliminary results on morphometric and allozyme variation in Belgian and French populations. *Nota Lepid. 1994;suppl.5:3–12.*

109. Dincă V, Runquist M, Nilsson M, Vila R. Dispersal, fragmentation, and isolation shape the phylogeography of the European lineages of *Polyommatus* (*Agrodiaetus*) *ripartii* (Lepidoptera: Lycaenidae). *Biol J Linn Soc Lond.* 2013; 109:817–829.

110. Horreo JL, Pelaez ML, Suárez T, Breedveld MC, Benoit H, Yann S-G, et al. Phylogeography, evolutionary history and effects of glaciations in a species (*Zootoca vivipara*) inhabiting multiple biogeographic regions. *J Biogeogr.* 2018; 45:1616–1627.

111. Trettin J, Agrawal S, Heinze J. Phylogeography of social polymorphism in a boreo-montane ant. BMC Evol Biol. 2016; 16:137. https://doi.org/10.1186/s12862-016-0711-3 PMID: 27334379

112. Vences M, Hauswaldt JS, Steinfartz S, Rupp O, Goesmann A, Künzel S, et al. Radically different phylogeographies and patterns of genetic variation in two European brown frogs, genus *Rana*. Mol Phylogenet Evol. 2013; 68(3):657–670. https://doi.org/10.1016/j.ympev.2013.04.014 PMID: 23632031

113. Todisco V, Gratton P, Zakharov EV, Wheat CW, Sbordoni V, Sperling FAH. Mitochondrial phylogeography of the Holarctic *Parnassius phoebus* complex supports a recent refugial model for alpine butterflies. *J Biogeogr.* 2012; 39:1058–1072.

114. Hegna RH, Galarza JA, Mappes J. Global phylogeography and geographical variation in warning coloration of the wood tiger moth (*Parasemia plantaginis*). *J Biogeogr* 2015; 42:1469–1481.

115. Kodandaramaiah U, Weingartner E, Janz N, Dalén L, Nylin S. Population structure in relation to host-plant ecology and *Wolbachia* infestation in the comma butterfly. *J Evol Biol.* 2011; 24:2173–2185. https://doi.org/10.1111/j.1420-9101.2011.02352.x PMID: 21745252

116. Kodandaramaiah U, Konvicka M, Tammaru T, Wahlberg N, Gotthard K. Phylogeography of the threatened butterfly, the woodland brown *Lopinga achine* (Nymphalidae: Satyrinae): implications for conservation. *J Insect Conserv.* 2012; 16:305–313.

117. Vandewoestijne S, Baguette M, Brakefield PM, Saccheri IJ. Phylogeography of *Aglais urticae* (Lepidoptera) based on DNA sequences of the mitochondrial COI gene and control region. Mol Phylogenet Evol. 2004; 31:630–646. https://doi.org/10.1016/j.ympev.2003.09.007 PMID: 15062799

118. Konvicka M, Zimmermann K, Klimova M, Hula V, Fric Z. Inverse link between density and dispersal distance in butterflies: field evidence from six co-occurring species. *Popul Ecol.* 2012; 54:91–101.

119. Mennechez G, Schtickzelle N, Baguette M. Metapopulation dynamics of the bog fritillary butterfly: comparison of demographic parameters and dispersal between a continuous and a highly fragmented landscape. *Lands Ecol.* 2003; 18:279–291.

120. Habel JC, Schmitt T. The burden of genetic diversity. *Biol Conserv.* 2012; 147:270–274.