Multiple Glacial Refugia of the Low-Dispersal Ground Beetle Carabus irregularis: Molecular Data Support Predictions of Species Distribution Models

Classical glacial refugia such as the southern European peninsulas were important for species survival during glacial periods and acted as sources of post-glacial colonisation processes. Only recently, some studies have provided evidence for glacial refugia north of the southern European peninsulas. In the present study, we combined species distribution models (SDMs) with phylogeographic analyses (using mitochondrial DNA = mtDNA) to investigate if the cold-adapted, stenotopic and flightless ground beetle species, Carabus irregularis, survived the Last Glacial Maximum (LGM) in classical and/or other refugia. SDMs (for both a western European and for a Carpathian subgroup) were calculated with MAXENT on the basis of 645 species records to predict current and past distribution patterns. Two mtDNA loci (CO1 and ND5, concatenated sequence length: 1785 bp) were analyzed from 91 C. irregularis specimens to reconstruct the phylogeography of Central and eastern European populations and to estimate divergence times of the given lineages. Strong intra-specific genetic differentiation (inter-clade ΦST values ranged from 0.92 to 0.99) implied long-term isolation of major clades and subsclades. The high divergence between the nominate subspecies and the Carpathian subspecies C. i. montandoni points to two independent species rather than subspecies (K-2P distance 0.042 ± 0.004; supposed divergence of the maternal lineages dated back 1.6 to 2.5 million years BP) differing not only morphologically but also genetically and ecologically from each other. The SDMs also inferred classical as well as other refugia for C. irregularis, especially north of the Alps, in southeastern Europe and in the Carpathians. The coincidences between the results of both methods confirm the assumption of multiple glacial refugia for the studied species and the usefulness of combining methodological approaches for the understanding of the history of low-dispersal insect species.


Introduction
The glaciers and permafrost grounds of the Quaternary ice ages shaped distribution patterns in biodiversity, as they caused species extinctions as well as the retraction of species' distribution ranges. Glacial refugia were important for species survival in glacial and interglacial periods and hosted many species which had recently been widespread across Europe [1]. Temperate and warmadapted species in particular are assumed to have colonised central and northern parts of Europe after the Last Glacial Maximum (LGM), originating from classical refugia such as the southern European peninsulas (Iberian, Apennine and Balkan) [1,2]. Cold-adapted species that had recently inhabited alpine and arctic habitats are assumed to have survived in the margins of southern European mountain chains [3]. In recent decades, palaeontological, palynological and phylogeographic studies have provided evidence for glacial refugia far north of the southern European peninsulas and the Alps (so-called cryptic refugia) [4,5].
In contrast to the species which were widespread across Europe [6,7], species with restricted and disjunct ranges in the European low and high mountain chains are poorly studied with regard to spatial patterns in genetic variation. However, some mountain plant species inhabiting the Alps and geographically relatively close mountain systems (e.g. the Central European mountain ranges, the Carpathians) show high genetic differentiation [8,9]. Zoological studies have covered arthropod species with high dispersal abilities (e.g. due to flight activity or ballooning): Butterflies [10] and spiders [11] show low level differentiations and indicate young and even post-glacial splits between populations from the Alps and adjoining high and low mountain ranges [12]. In contrast, low-dispersal species are likely to show stronger genetic differentiation patterns due to the absence of gene flow (e.g. flightless carabids; [13,14]); however, the population histories of these species have been studied only insufficiently to date.
New study methods such as species distribution models (SDMs) have become useful tools to predict potential habitat distribution of species under current climate conditions and -when analysed in conjunction with palaeoclimatic data -to project species' past distributions [15]. Many studies detected classical as well as cryptic refugia using past predictive models (e.g. [16,17]). Unaffected by incomplete sampling, the extinction of genetic variants and largescale range shifts of species, past predictive models can contribute to the localization and sizing of species' glacial refugia [18]. Consequently, past predictive models and phylogeographic analyses can be used complementarily to develop reasonable scenarios of species' range retractions during glaciations and of range expansions during interglacial stages.
In this study we applied a combined approach using modelling and of phylogeographic analysis of two mitochondrial DNA loci (CO1 and ND5) to investigate the history of the flightless ground beetle species, which occurs in a disjunct range across high-and low-altitude mountains of Central and eastern Europe: Carabus irregularis. We aimed to test (i) whether the cold-adapted species survived the LGM in classical and/or cryptic refugia and (ii) whether the low-dispersal species shows old phylogenetic groups and high intra-specific variation. Since the CO1 locus has been successfully used in barcoding [19], we also checked (iii) whether the phylogeny of Carabus irregularis is consistent with the intraspecific taxonomy of the species.
We addressed the following specific study questions: (1) Where did C. irregularis find suitable climate conditions (refugia) to survive the LGM? (2) Does the species' genetic differentiation pattern provide molecular evidence for low dispersal and for the detected potential glacial refugia? (3) Is the phylogeny of C. irregularis congruent with its taxonomic subspecies delineation?

Study species
Carabus (Platycarabus) irregularis Fabricius, 1792 occurs in cold and mesic habitats in montane to subalpine altitudes (e.g. beech forests and high mountain meadows). Suitable habitats are predominantly on limestone and are rich in snails, the preferred food of C. irregularis. The disjunct distribution of C. irregularis ranges from the Teutoburg Forest in the north to the Dinaric Mountains in the south and from the Carpathians in the east to the Jura Mountains in the west ( Fig. 1) [20]. In the taxonomic literature, three subspecies are well-accepted even though different names are sometimes used [21]. In this study we follow Casale and Kryzhanovskij [22] and use the following binomials: C. i. irregularis s.str. in Central Europe including the northern Alps, C. i. bucephalus Kraatz, 1879 in the southeastern Alps and the northern Balkan Peninsula, and C. i. montandoni Buysson, 1882 in the Carpathians (Fig. 1). While the distribution ranges of C. i. irregularis and C. i. bucephalus pass into each other in the southeastern Alps, C. montandoni is geographically separated from the other two subspecies (see subspecies affiliation in Table S1). The species is declining in some regions [23,24] and already extinct in some parts of its former distribution range (e.g. in Belgium [25]).

Species distribution modelling
In order to predict the potential habitat distribution of the study species under current and past climate conditions, we compiled 1,005 data points of observed occurrences of C. irregularis from the literature and databases (e.g. [26][27][28][29][30][31][32][33]), from museum and private entomological collections as well as our own field surveys (see complete list of records in Table S1). Since spatial clumping of species records due to different sampling intensities may bias SDMs, redundant localities were excluded by the model software, and the final locality input included 645 records. In addition to the geographic separation, there are also differences in the altitudinal distribution pattern of the two western European subspecies (C. i. irregularis and C. i. bucephalus) and the Carpathian subspecies (C. i. montandoni). According to the compiled species records, C. i. irregularis and C. i. bucephalus occur in high as well as low altitudes in the western part of its distribution range and on average in lower altitudes (altitudes as low as 250 m a.s.l.; mean altitude of western records = 774 m a.s.l.; standard error (SE) = 13.71; confidence interval (CI) = 26.87; see Fig. S1). From the eastern part of the range (the Carpathian mountain system), however, the species was scarcely recorded from altitudes below 600 m a.s.l. (mean altitude of Carpathian records = 912 m a.s.l.; SE = 45.50; CI = 89.19; see Fig. S1). Due to geographically separated distribution ranges and non-overlapping confidence intervals of the altitudinal distribution of the subspecies, we could not exclude different habitat preferences of the Carpathian subspecies. Thus, species records were split into two sets before modelling: one included occurrences of the Carpathian subspecies, the other comprised occurrences of the two other subspecies from the western part of the distribution range (Central Europe and Balkan). Climate data (monthly temperature and rainfall values) was generated according to Hijmans et al. [34] to 19 bioclimatic variables (BIO1-BIO19; Table S2) which representing factors of climate which are particularly significant in determining species' distributions [18]. For carabids, temperature and humidity variables have been demonstrated to be the most important environmental factors influencing habitat selection [35].
Current climate data (i.e. period 1950-2000) was downloaded with a spatial resolution of 30 arc s (,161 km) from the WorldClim database (http://www.worldclim.org/). Past (Last Glacial Maximum = LGM; 21,000 BP) climate data downscaled from the MIROC general circulation models to the required spatial resolution (2.5 arc min) was provided by the Paleoclimate Modelling Intercomparison Project Phase II (http://pmip2.lsce. ipsl.fr/). Since our SDMs were to include only predictors with high explanatory power and ecological relevance [36,37], we included only the variables which are most likely to influence the occurrence of C. irregularis. Therefore, we extracted all 19 bioclimatic variables at coordinates of species records and performed a principal component analysis (PCA). To reduce multicollinearity, we compiled a correlation matrix and excluded all strongly correlated variables (Pearson's r 2 .0.75; [38]). Final models included five uncorrelated predictor variables: 'Isothermality' (BIO3), 'temperature seasonality' (BIO4), 'minimum temperature of coldest month' (BIO6), 'precipitation of wettest month' (BIO13) and 'precipitation of driest month' (BIO14) -climatic measures reflecting habitat conditions which are likely to affect the occurrence of C. irregularis [30,31]. The bioclimatic layers were cut to cover the recently documented distribution range of the species sufficiently but not too generous, since large layer extents could have a negative effect on SDM results [39]. The layers used ranged from 30 to 60 decimal degrees North and from 0 to 30 decimal degrees East.
We used the maximum entropy approach [40], which repeatedly performed better than other methods using species presence data in comparative studies [41,42]. Three sets of SDMs -one with the western record set (hereafter western distribution model), one with only Carpathian records (hereafter Carpathian distribution model) and one model including all records together (hereafter entire distribution model) were generated with MAX-ENT 3.3.3 (http://www.cs.princeton.edu/schapire/maxent), with automatically sampled random background points and species records split into 75% model training and 25% model evaluation data. Model quality was evaluated by a bootstrap validation and the Area Under Curve method (AUC.0.7; [43]). The logistic output of 100 model replicates was transformed to species maps showing the potential distribution of C. irregularis under current climate conditions and under LGM conditions. Species maps were processed in DIVA-GIS v7.1.7 (http://www.diva-gis.org). We compared the Maxent outputs of the three model sets using the software ENMtools [44]. Niche overlap was quantified by two metrics of niche similarity introduced by Warren et al. (2008) [45]: D (Schoener's statistic for niche overlap; [46]) and I (derived from Hellinger distance). Both measures range from 0 (no niche overlap) to 1 (identical niches). In addition, a niche identity test was run to evaluate niche similarity between the western and the Carpathian distribution model [44].

Sample collection and molecular methods
In total 91 specimens of the study species were collected from 23 localities across its entire distribution range (Table 1, Fig. 1). In spite of conscientious efforts and communications with numerous coleopterologists, no fresh beetles were available from certain parts of the distribution range (possibly due to the recent decline mentioned above). Adult beetles were collected from hibernation in autumn or caught by pitfall trapping during the activity periods. The specimens were preserved at -80uC or in pure ethanol at -20uC until DNA extraction.
Mitochondrial genomes within the genus Carabus can be influenced by horizontal gene flow (secondary contact after speciation) [47][48][49]. In natural habitats, C. irregularis seems to produce hybrids with related species (members of the same subgenus Platycarabus: C. creutzeri and C. depressus) [50]. Therefore, we used six individuals of the latter species from three different sites as outgroups for rooting phylogenetic trees and for checking the reliability of mitochondrial DNA (mtDNA) as a phylogenetic marker. In comparative studies, the two different mitochondrial genes studied, CO1 and ND5, have outperformed nuclear markers in terms of giving a reliable estimates in timing of splits and phylogenetic reconstruction because of their maternal heritability [51]. Moreover, CO1 is a valuable barcoding marker in animals [19] which makes it possible to study the performance of such a marker system on an infraspecific level.
Genomic DNA was extracted from one femur of each specimen by CTAB (Cetyldimethylethyl-Ammoniumbromid) lysis [52], where we reached a DNA concentration of min. 200 ng/ ml. Cycling conditions for CO1 started with an initial denaturation for 5 min at 94 uC, followed by 35 cycles of denaturation for 45 at 94uC s, annealing for 45 s at 50uC and extension for 1 min at 72uC. Final elongation was performed for 8 min at 72uC. For ND5 the following PCR program was used: initial denaturation at 94 uC for 5 min, 35 cycles (1 min at 94u, 1 min. at 50uC, 2 min at 70uC), then final elongation for 7 min at 70uC.
PCR products were visualised on 2% agarose gels, dyed with Roti-Safe Gelstain (Carl Roth GmbH, Karlsruhe, Germany) and then purified with the GeneJET PCR Purification Kit (Fermentas GmbH, St. Leon-Rot, Germany). The manufacturer's protocol was slightly modified: only 50 ml elution buffer were used for final elution to ensure a minimum DNA concentration of 100 ng/ ml. Purified PCR products were again checked on agarose gels. Sequencing (both forward and reverse strands) of PCR products was carried out at the laboratory of LGC Genomics (Berlin, Germany) using an Automatic Sequencer 3730 xl (Applied Biosystems, Foster City, USA).

Sequence alignment and phylogenetic analysis
Forward and reverse sequences were assembled and manually corrected using the software Geneious v4.8.5 [55]. Sequences were aligned and checked for reading frame errors in MEGA5 [56]. Prior to phylogenetic analysis, MrModeltest v2.3 [57] was used to identify the best nucleotide substitution model for each mtDNA fragment (CO1 and ND5) ranked by Akaike Information Criterion (AIC). For each gene, we performed an independent run of 60610 6 generations, sampling one tree every 1,000 generations. Based on the selected evolutionary models, we investigated phylogenetic relationships within the species C. irregularis (including two outgroups) using Bayesian inference (BI), maximum likelihood (ML) and maximum parsimony (MP). The BI tree was calculated by MrBayes v3.2.1 [58] running 2610 6 generations, while reaching an average standard deviation of split frequencies,0.01.
The ML method was conducted in RAxML v7.3.0 [59] performing 10,000 thorough bootstrap replicates. MP analysis was carried out in MEGA5 [56] using 5,000 bootstrap replicates. Divergence times were estimated by using strict molecular clocks in BEAST v1.7.1 [60] considering gene-specific mtDNA substitution rates for the two studied loci (1.45% My 21 for CO1, and 1.59% My 21 for ND5) based on rates determined for other carabid beetle species [51,61]. Population size through time was applied by a Yule speciation model, and Markov-Chain-Monte-Carlo (MCMC) chain length was set to 200610 6 , whereas other parameters remained in default settings.

Intra-specific divergence and genetic diversity
Measures of DNA sequence variation within different hierarchical levels (populations, subclades and clades) and in all sequences together (total) were computed for concatenated sequences and for CO1 and ND5 separately using DnaSP v5.10.01 [62]. Variation measures included the number of polymorphic/segregating sites (S), number of haplotypes (H), haplotype diversity (h) and nucleotide diversity (p). In addition, we performed an analysis of molecular variance (AMOVA) using ARLEQUIN v3.5 [63] to estimate the level of genetic differentiation (W ST ) among populations, subclades and clades. We used the median-joining approach [64] implemented in the NET- WORK Software v4.6.1.0 (http://www.fluxus-engineering.com) to illustrate phylogenetic and geographic patterns in haplotype diversity. Intra-specific genetic divergence in CO1 -both within and between major clades -and inter-specific genetic divergence in the studied outgroup species of the subgenus Platycarabus were calculated by pairwise comparisons (4,005 intra-specific pairs and 554 inter-specific pairs) using Kimura's two parameter (K-2P) model in MEGA. Teutoburg Forest, Weser Mountains and Harz) and also in the northern Carpathians (Tatra Mountains). Low suitability (0.3,p,0.5) was attributed to the eastern Carpathians (Romania, Ukraine), the Ardennes (France, Belgium, Luxembourg), large parts of Germany and also up to Denmark. Areas with high suitability rates were always surrounded by areas with lower rates, resulting in a disjunct distribution of suitable habitats. Thus, the model was widely consistent with the documented actual distribution of C. irregularis. However, the study species is not known from the regions with low suitability rates, with the exception of the Romanian and Ukrainian Carpathians. The Carpathian distribution model for current climate conditions (Fig. 2b) displayed the most suitable conditions in parts of the Alps and several high-altitude regions in the Carpathian mountain system. The past predictive model of C. irregularis using the western European data set showed the most suitable conditions (p.0.7) during the LGM along the edges of the large Alpine glacier: in the southwestern Alps, the Vosges and Jura and in the Ligurian part of the Apennines. Well suited climatic conditions were also shown in a large area on the Balkan Peninsula (large parts of the Dinaric Alps) and in a small part of the northwestern German low mountain range. Moderately suitable regions were found in the central Apennines, in the eastern parts of the Alps and in the western Balkan. In the Carpathians, models for LGM conditions showed only small regions with low suitability rates (see arrows in Fig. 2c).

Current and past predicted distribution
The past distribution model for the Carpathian data set (Fig. 2d) illustrated large areas with highly suitable climate conditions in the southern parts of the Eastern Alps and also high suitability rates for the Tatra Mountains and for some parts of the Southern and the Serbian Carpathians.

Phylogenetic analyses and divergence time estimation
For both loci, the best evolutionary model identified by the AIC was the general time reversible model with Gamma distributed rate variation across sites and a proportion of invariable sites (GTR+G+I) [65]. Since there has been concern about estimation problems in this model type due to interactions between the proportion of invariable sites and the Gamma distribution [59], we also checked the model without +I and chose this simpler model for the ML analysis in RAxML.
BI, MP and ML yielded highly congruent phylogenetic trees with the same major nodes and branching order. The majorityrule consensus tree (Fig. 3) displays the Bayesian posterior probability (BPP), parsimony bootstrap percentages (PB) and maximum likelihood bootstrap percentages (MLB), where BBP $ 0.5 and bootstrap values (PB and MLB)$50% represent wellsupported nodes [66]. In addition, nodal ages and 95% confidence intervals are represented for each node.
In the phylogenetic tree, all basal nodes had very high support (BPP$0.98; PB$93; MLB$65) and branches represent two major  Fig. 3). Assumed divergence times spanning the last 2 million years (Fig. 3) indicated that major clades and subclades within C. irregularis were formed long before the last glacial period. The separation of C. irregularis from other species of the subgenus Platycarabus was found  Table 1 for population abbreviations and Figure 1 for clade abbreviations and colour codes. doi:10.1371/journal.pone.0061185.g003 to have taken place earlier than the separation of the species used as outgroups (C. creutzeri and C. depressus). Additionally, the split into the two major lineages within C. irregularis appears to be older than the split of the two species C. creutzeri and C. depressus. Only one of the three subspecies was monophyletic (C. i. montandoni), while C. i. irregularis and C. i. bucephalus are paraphyletic.

Intra-specific diversity and geographic distribution of genetic diversity
The overall haplotype diversity (h = 0.882 6 0.0006) and nucleotide diversity (p = 0.01760.002) were high, with clade A showing higher haplotype and nucleotide diversity (h = 0.901; p = 0.016) than clade B (h = 0.838; p = 0.009; Table 2). In general, haplotype and nucleotide diversity tended to be higher in populations from the southern part of the distribution range (Fig. S3). Within clades, subclade As showed higher diversity than Ae and the subclades Bsl and Bs are more diverse than subclade Bn (see Table 3 for a summary of computed diversity measures of different hierarchical levels). The mean intra-specific variation in CO1 (K-2P; d = 0.027 6 0.004) was also very high. The interclade K-2P divergence was 0.042 6 0.004 and the frequency distribution of intra-specific genetic divergence in CO1 within and between major clades showed a clear pattern (Fig. 4). On the population level, K-2P distance values ranged from 0.0 to 0.016 between populations within clade B to 0.022 between populations within clade A. Between populations of clade A and B, genetic distance ranged from 0.040 to 0.044 (Table S5).
In addition, we detected a significant differentiation structure among all populations (W ST = 0.935, p,0.001) and no shared haplotypes between the clades, subclades and subspecies but between populations within clades and subclades. Additionally, AMOVAs for different hierarchical groups (subclades, clades, subspecies) revealed significant genetic variance among these hierarchical groups, among populations and also within populations of the hierarchical groups. Most of the diversity was observed among subclades (88.1%), among clades (72.3%) and among subspecies (63.7%), while low percentage of variance was detected within populations (3.1 to 4.8%; Table 3). Differentiation between the single populations of clade B ranged from lower values in geographically closer populations (W ST = 0.10) to higher values in geographically distant populations (W ST = 0.99). Populations of clade A also showed very high W ST values (.0.9; Table S5), even though the sampled populations were geographically closer than those from clade B. Inter-clade W ST values comparing populations of clade A and clade B ranged from 0.92 to 0.99. The haplotype network (Fig. 5) supports these findings and illustrates a similar pattern: high genetic distance between the two geographically distant (,1,500 km) major clades, but also large genetic distance and a high number of substitutions between the geographically closer (,300 km) Carphathian subclades (Ae and As). In contrast, populations within subclade Bs showed very low genetic distance, although geographic distance is quite large (between Bla/Sw I-III/ Jur I-II and SloI ,900 km).

Potential glacial refugia
Our past predictive models showed several independent mountainous regions across Central and southern Europe with highly suitable climate conditions for C. irregularis which could have acted as glacial refugia during the LGM. Some of these potential refugia (e.g. in the Dinaric Alps) conformed to classical refugia on the southern European peninsulas [1,67,68] and other refugia at the margins of the Alps and in neighbouring mountain ranges (e.g. the Jura, the Black Forest), some of which have quite recently been recognised as refugia [69][70][71]. In addition, our model pointed to refugial areas far north of the southern European peninsulas and the Alps. So-called cryptic refugia have been assumed for a growing group of other species (plants as well as animals) [5,72].
While the southern European peninsulas are typical glacial refugia of current lowland plant and animal species [2,73], mountain species potentially persisted in refugia along the border of the Alps when higher altitudes were covered with ice sheets Figure 4. Frequency distribution of intra-specific genetic divergence in CO1 -within (light grey) and between major clades (dark grey). Pairwise distances (4005 intra-specific comparisons within C. irregularis) were calculated using Kimura's two parameter (K-2P) model. doi:10.1371/journal.pone.0061185.g004 [69,74]. According to our results, C. irregularis seems to have survived in refugia of both low-and high-altitude species. Our models projected relatively large and consistent areas of suitable and moderately suitable conditions during the Last Glacial Maximum (comprising the maximum extension of the continental glaciers), suggesting that C. irregularis was relatively widespread during glaciations. In general, numerous mountain species (even though mostly defined as species with their main distribution above the timber line) are assumed to have survived in several spatially restricted and isolated refugia rather than in single extensive and consistent refuge areas ( [75] and references therein; [76]). Consequently, C. irregularis, a cold-adapted montane to Table 2. Diversity statistics based on concatinated CO1+ND5 sequences (1785 bp) of C. irregularis populations (sampled with more than one specimen) and in subclades and clades (as determined by the phylogenetic tree).  subalpine species, can be expected to have survived the last ice age in multiple spatially restricted refugia within the projected potential distribution area in the Alps and some neighbouring mountain ranges. It appears less important whether there is a significant difference between the historical and the recent distribution range. The fact that the study species tended to have cryptic refugia north of the Alps seems to be of greater importance. While cold-tolerant species, e.g. some mammalian herbivores, have long been believed to have survived in Central European steppe tundra surrounding the glacial ice-sheets, cryptic refugia were surprisingly also found for more temperate animal and tree species [72,77]. Topographically sheltered sites are assumed to have provided suitable climate conditions and stable microclimates for species with differing habitat requirements [3,72]. Species distribution modelling is a relatively new but already frequently used method for inferring species potential distribution ranges from occurrence data. The method has been criticised for including abiotic factors only and not taking biotic parameters such as competition and dispersal into account, which might explain discrepancies between the documented range and the projected species range [15,16]. Nonetheless, the high model quality of our simulations (AUC = 0.92 and 0.98) represented high consistence between actual and predicted occurrence. Another uncertainty, especially of past predictive modelling, is the accuracy of past climate data. We used data processed by the downscaling method, which is currently the best large-scale data available but still requires further validation [18]. Additionally, the modelling approach is based on niche conservatism and does not consider changes in habitat requirements. Nevertheless, past predictive models are useful tools to visualise species distribution patterns and potential glacial refugia. More traditional approaches to identifying refugia are pollen and (sub-)fossil analyses (e.g. [78]). Since no subfossils are known for Platycarabus species (cf. [79]), we chose a molecular method (mtDNA analysis) to supplement our findings from modelling and to review assumptions concerning potential glacial refugia.

Genetic evidence for glacial refugia
Phylogenetic analyses (using BI, MP and ML approaches) revealed two well-supported and spatially congruent major clades: one Central European (B) with three subclades and one Carpathian clade (A) with two subclades. The ages of all major phylogenetic splits predating the last ice age, the absence of shared haplotypes and high genetic intra-specific differentiation between clades as well as subclades implied that the ancestors of phylogenetic (sub)clades in C. irregularis survived many glacialinterglacial cycles isolated from each other. The Central European clade consisted of three independent and genetically distinct Central European subclades: one subclade with a more southern distribution (Bs: with populations from southern Germany, France and Slovenia), one Slovenian subclade (Bsl) and another one including haplotypes from the northern part of the species' distribution range (Bn: central to northern Germany). Thus, the assumption of multiple refugia of C. irregularis in Central Europe inferred from past predictive models was supported by our phylogenetic results, with at least two refugia located close to the Alps or on the Balkan Peninsula and one northern refugium, potentially in central Germany. However, phylogenetic analyses indicated the existence of independent Carpathian refugia, which were hardly evident from the western European distribution model for LGM conditions (Fig. 2c, arrows). The past predictive model using the eastern records showed several potential refugia in the Carpathian mountain system (Fig 2d). These potential refugia coincide with the finding of two subclades of clade A one from the eastern (Ae) and one from the southern Carpathians (As), with divergence times dating back more than one million years BP.

substitutions
As well as refugia in northern Central Europe, Carpathian refugia are not accepted as classical refugial areas and phylogeographic studies including the Carpathian region are still rare. However, since the Carpathians have faced a different glacial history with more local glaciers and larger areas of suitable alpine habitats than the well-studied Alps [9,80], investigation of the southeastern European mountain range appears likely to be very interesting. Some recent studies (on mountainous plant and animal species) detected low-level differentiation between the Carpathians and the Alps and assumed the Carpathians to have been colonised by adjacent refugia e.g. in the western Alps [8] or the Balkans [71]. In contrast, an increasing number of phylogeographic studies have found very old and divergent phylogenetic lineages, suggesting a major refugium for plants [9,80] as well as for vertebrates [81][82][83][84] in the Carpathian region. In congruence with our phylogeographic results for C. irregularis, alpine plants, namely Hypochaeris uniflor and Campanula alpina, showed distinct and strongly differentiated genetic groups within the Carpathians, also pointing to long-term isolation and restricted gene flow between several areas within the mountain system [9,80,85]. While Ronikier et al. [80] discovered the strongest differentiation between western and eastern Carpathian samples, results of Mráz et al. (2007) [85] resembled our findings in that they also showed high differentiation between southern and eastern Carpathian populations. Examples from animals also revealed genetically distinct and longterm isolated phylogenetic groups in the Carpathians [81][82][83][84].

DNA barcoding and systematic delineation
The phylogenetic results of our study agreed only in part with subspecies systematics: Whereas the Carpathian subspecies C. i. montandoni was monophyletic (clade A), the other two subspecies (C. i. irregularis and C. i. bucephalus) were paraphyletic (clade B) and we suggest to rank the populations belonging to clade B as one subspecies. The high proportion of polymorphic sites within ND5 (7.70%) and also within the CO1 loci (8.52%), which is commonly used as a barcoding marker [19] -also for ground beetles [86], questions the conventional taxonomy. DNA barcoding using the CO1 locus for rapid species delineation [87], implies that interclade K-2P distances higher than 2.8 to 3.4% indicate species rather than subspecies or populations. Following these thresholds, the divergence of 0.042 6 0.004 between clade A (taxonomic subspecies C. i. montandoni) and B (irregularis s. str.) points to independent species rather than subspecies of the same species. Although the commonly used Kimura's two-parameter substitution model (K-2P) has been repeatedly criticised [88,89], estimates using other substitution models showed similar results for species delineation [89]. For C. irregularis, the calculation of other (simpler as well as more complex) genetic distance metrics shows a consistent pattern with K-2P (p distance between clade A and B = 0.042; Tamura-Nei distance between clade A and B = 0.044).
However, due to the fact that substitution rates can vary even within the same species [90] it is necessary to exercise caution when defining a species from nucleotide rates alone.
Using morphological and molecular data Casale et al. [91] reconstructed the phylogeny of the subgenus Platycarabus. Following their phylogeny C. irregularis shows the strongest differentiation from all other species within the subgenus (basal split) and C. depressus and C. creutzeri are the second major diverged species pair. Surprisingly, the split between the two main lineages of C. irregularis we found is older than the divergence between the latter species pair. This underlines both (i) the possible species status of the Carpathian populations and the (ii) need to incorporate a population-based approach to understand the phylogeny of these highly differentiated ground beetles with their ancient lineages, even below or at the species level.

Conclusion
Past predictive modelling and phylogenetic analyses acted as supplements and imply that C. irregularis survived the last glacial periods in long-term isolated classical refugia on the edges of the Alps as well as in other refugia (in Central Europe and some parts of the Carpathian Mountains). Altogether our results indicate that the Carpathian subspecies of C. irregularis differs not only in terms of its geographical distribution but also genetically and (due to differing habitat preferences as revealed from the poorly overlapping SDMs) ecologically from rest of the species. Thus we conclude that C. irregularis comprises at least two evolutionarily significant units (ESUs, sensu Moritz, [92]) indicating that the species has a very interesting history, which should be investigated in more detail by further phylogeographic analyses -also in the framework of the Platycarabus group.   Table S4 Summary of statistical analyses on measures of DNA sequence variation within clades, subclades and in all sequences together (total) for CO1 and ND5 separately: N = sample size, S = number of polymorphic/segregating sites, Eta = total number of mutations, k = average pairwise nucleotide difference per sequence, H = number of haplotypes, h = haplotype diversity, p = nucleotide diversity, hG = mutation parameter per sequence).

(XLSX)
Table S5 Genetic differentiation (W ST ; in the upper part of the matrix) and intra-specific genetic divergence (K-2P distances; in the lower part of the matrix) of Carabus irregularis. K-2P distances greater than the threshold value assumed for species delineation by Wiemers and Fiedler [87] are printed in bold. (XLSX)