The specimens used for this work were mainly collected in the Pyrenean and Cantabrian mountain chains (Fig. 2). Further specimens included were from the Massif Central (France), Italy, Bosnia Herzegovina and Turkey. The respective collection sites are given in Table 2.

Directly after collection in field, the specimens were transferred to 2 ml plastic microtubes with sealed screw caps, filled with 95% ethanol to preserve the specimens. The tubes were later stored at –20° C. The specimens of each collection site were sorted by morphospecies. One individual of each morphospecies was used for DNA analysis. A code was then given to each of these specimens (e.g. ‘Cr2’, Table 2). These codes are referenced throughout this work. Spare specimens remained stored as described above. This procedure was carried out for all collection sites. In total, 26 specimens were sampled and processed for DNA analysis. The specimens are deposited at the Stuttgart State Museum of Natural History (SMNS).

DNA extraction and purification were carried out with the DNeasy Blood & Tissue Kit (50) (QIAGEN, Hilden, Germany) following the manufacturer’s instructions. The DNA extraction was non-destructive and processed with whole specimens. The areas between the head and the thorax and between the thorax and the abdomen were slightly opened to allow better digestion by the proteinase K. Overnight sample incubation for DNA extraction was performed with the Heating ThermoMixer MHL 23 (Ditabis, Pforzheim, Germany). Subsequently, the DNA concentration of each sample was measured using the NanoPhotometer® N60 (IMPLEN, München, Germany) to confirm a successful extraction.

Four gene fragments were targeted for sequencing including two mitochondrial and two nuclear genes: “cox1” cytochrome c oxidase subunit 1 – mitochondrial (CO1); “rrnl + tRNA-Leu + nad1” 5’ end of the large ribosomal 16S unit + tRNA-Leucine gene + 3’ end of the NADH dehydrogenase subunit 1 – mitochondrial (16S); “LSU” large ribosomal subunit – nuclear (28S); “SSU” small ribosomal subunit – nuclear (18S). The primers used are listed in Table 1. For some samples, the amplification of the whole fragment of CO1 and LSU fragments failed. In those cases, primer pairs targeting shorter but partially overlapping fragments were used. DNA amplification was carried out with the Taq PCR Master Mix Kit (250) (QIAGEN, Hilden, Germany) using the Labcycler Basic (SensoQuest, Göttingen, Germany). A standard program was first used for all PCRs. In case of an unsuccessful run (no–, insufficient– or ambiguous results shown in the control gel), the cycling programs were slightly altered to achieve better amplification results.

The PCR products were then controlled with a gel electrophoresis. Therefore a 1%-agarose gel (1:100 agarose / TAE buffer 1×) with GelRed® (Fremont, CA, USA) was run at 100 V in a Mupid® One Electrophoresis System, Advance (Mupid CO. LTD., Tokyo, Japan) for 25 min. The gels were then photographed under UV light using the Pentax TV Zoom lens 8–48 mm 1:1.0 (Ricoh Co. Ltd. Operations, Tokyo, Japan) and the BioDocAnalyze-Software (Analytik Jena GmbH, Jena, Germany). PCR product purification was conducted with the QIAquick PCR Purification Kit (250) (QIAGEN, Hilden, Germany). In preparation for sequencing, 5 µl of the purified PCR product and 5 µl of the respective primer were added to a 5 ml centrifuge tube. The same primer aliquots were used to reduce the possibility of contamination after the PCR product was controlled with a gel. The samples were then sent to the Macrogen laboratory Europe B.V. (Amsterdam, Netherlands) for sequencing.

The raw sequences were processed with GENEIOUS PRIME® 2020.2.2 (https://www.geneious.com). The sequences were cleaned and aligned with „Geneious Alignment“ (global alignment with free end gaps, cost matrix: 65% similarity) and primer sequences were trimmed using the „trim primer“ function. Consensus sequences were aligned with MUSCLE v.3.8.425 (R. C. Edgar, www.drive5.com/muscle/). Additional sequences of Pterostichus specimens available at Genbank (www.ncbi.nlm.nih.gov/Genbank) were added to the alignments. All specimens included in this work are listed in Table 2. The alignments were then exported to BIOEDIT v.7.2.5 (Hall 1999). Single sequences within the alignment were brought to the same length. Therefore, parts of the 5’ or 3’ ends of some sequences were cut or filled with a placeholder character which was then recognized as missing information by the following programs used in the phylogenetic analysis.

The cleaned alignments were exported in NEXUS format. For the combined analysis, the alignments were first assembled in a data matrix using SEQUENCE MATRIX v.1.8 (Vaidya et al. 2011). Tree reconstruction was performed by maximum likelihood analysis for all single genes as well as the combined matrix. W-IQ-TREE 1.6.12 (Trifinopoulos et al. 2016), available at the IQ-TREE web server http://iqtree.cibiv.univie.ac.at, was used for maximum likelihood analysis including ultrafast bootstrap (Hoang et al. 2017). The MODELFINDER tool was applied to determine the best fitting substitution model beforehand (Kalyaanamoorthy et al. 2017). The chosen models were TIM + F + I + G4 for CO1, TPM2u + F + I + G4 for 28S, TIM2 + F + I for 16S, JC for 18S and GTR + F + I + G4 for the combined alignment (Jukes and Cantor 1969; Kimura 1981; Tavaré 1986; Posada 2003). In addition, a Bayesian inference was performed for the combined matrix using MRBAYES 3.2.7a (Ronquist and Huelsenbeck 2003) applying the Markov chain Monte Carlo algorithm (MCMC) altered by Geyer (1991). The output was visualized with FIGTREE v.1.4.4 (http://github.com/rambaut/figtree/) and edited using Adobe Illustrator v.26.0.3 (https://adobe.com/products/illustrator). For the combined analysis, a chimera was created for P. riparius (Dejean, 1828) (JF888281 + EU142445). The determination of the respective specimens was trusted. Only the topology of the combined analysis is shown in the results section (Fig. 3). For the single gene analyses, a summary is given in Table 3.

After DNA extraction the specimens were glued on rectangular cards for morphological study. For male individuals, the genitalia were removed beforehand and glued beside the specimen. Species determination was conducted regarding the currently valid species list of Cryobius published in the “Catalogue of Palearctic Coleoptera V1” (Bousquet 2017). The identification of the Pyrenean and Cantabrian species was mainly performed using the dichotomic key published in the “Faune de France Vol. 95” (Pupier 2011) and the “Faune de France Vol. 40” (Jeannel 1942) for male genital characteristics. The Turkish specimen was tentatively identified using species publications with type locality near the collection site, as there is no determination key available that includes all Palearctic Cryobius species. Additional information provided by the DNA analysis or distribution patterns of Cryobius species (Serrano 2013; Coulon and Pupier 2014; Bousquet 2017) was used to validate the morphological determination. When possible, determination was carried out to subspecies level. However, in some cases, the variability of external characters even within populations (Pupier 2011) or missing information concerning characters of female genitalia did not allow for an unambiguous classification at the subspecies level.

Species identification of non-Cryobius specimens was conducted with the “Käfer Mitteleuropas – Band 2, Adephaga 1” (Müller-Motzfeld 2004).

Our data support the monophyly of Cryobius with a bootstrap value (BV) of 100 and a Bayesian probability of 100 (BP) in the combined tree (Fig. 3). This monophyly is also recovered by the CO1 and the 16S single gene analyses (BV = 100 and 78, respectively; Table 3). The two conserved markers 28S and 18S do not show a distinct Cryobius clade. However, it must be said that testing the monophyly of Cryobius was not the scope of this work and the type species of Cryobius (C. ventricosus (Dejean, 1831)) was not included, although two species belonging to the ventricosus group sensu Ball (1966) (C. caribou Ball, C. riparius Dejean) were included. Also, a much more comprehensive sampling of Pterostichus subgenera and species would be required to achieve this aim.

No support was found for the synonymized subgenera Haptoderus and Pyreneorites. Neither the combined – nor the single gene phylogenies showed a distinct clade for either Haptoderus – or Pyreneorites species (Fig. 4, Table 3).

The Pyrenean and Cantabrian specimens form a monophyletic clade together with P. pumilio. The local Pyrenean species P. infimus is, regarding the combined phylogeny, grouped in one clade with P. pumilio (BV = 87, BP = 95). This arrangement is also found by 16S (BV = 70) and 18S, but here with low support. In the CO1 analysis, P. infimus is placed as sister to P. pumilio (BV = 79), which is not the case in the 28S phylogeny. According to the combined phylogeny, the remaining species with an exclusively Pyrenean or Cantabrian distribution are grouped in three separate clades.

The widely distributed species P. pumilio does not, according to the combined tree, form a monophyletic group since the Pyrenean P. infimus belongs to the clade. A monophyly for P. pumilio is only recovered in the CO1 topology (BV = 94). The clade of Pterostichus pumilio is further divided into two subclades. The first is a clade of specimens from the Massif Central (Cr20 + Cr27.2) and Germany (BV = 99) and BP = 100 in the combined tree). The clade is also recovered and well supported within the CO1 tree (BV = 92) (Table 3: “pumilio excl. Pyrenees/Cantabria”). No statement concerning the support for this clade is possible for the three other markers taken independently (16S, 28S, and 18S) as for those gene fragments only one sequence was available for this group. The second clade comprises Pyrenean and Cantabrian P. pumilio with a BV of 61 and a BP of 88 in the combined tree (Table 3: “Pyreneo-Cantabrian pumilio”). It has low support in the 28S and 16S trees (not tested in 18S). The clade was not recovered in the CO1 phylogeny.

Except for P. infimus, all Cryobius specimens with an exclusive Pyrenean distribution are arranged in two sister clades. The clade containing P. abaxoides and P. colasi is very well supported (BV = 97, BP = 100). The adjacent clade contains P. pusillus and P. amoenus with a relatively good support (BV = 75, BP = 82). Within those Pyrenean Cryobius, two subclades are well supported: one gathering all specimens of P. abaxoides (BV = 100, BP = 100) and one including all specimens of P. pusillus (BV = 100, BP = 100). Both are based on two specimens each.

Another well supported clade includes the P. cantabricus – and P. aralarensis groups (BV = 100, BP = 100), as well as two P. cf. subiasi (Ortuño et Zaballos, 1992) specimens (further discussed below). This clade will be referred to as the ‘Cantabrian clade’ (Table 3: “Cantabrian Cryobius ex pumilio”; including all Cantabrian specimens, P. pumilio excluded). The Cantabrian clade is also found in the CO1 tree (BV = 96), the 28S tree (BV = 74) and the 16S tree (BV = 81). It is not recovered in the 18S single phylogeny.

Two Spanish specimens (Cr9.2: female, Cr10: male) included in the Cantabrian clade were determined as P. (C.) cf. subiasi. Based on their outer morphology and their position within the tree, it was assumed that they are the same species. The uncertainty of the determination at species level was due to several ambiguous clues. The morphology of the aedeagus of specimen Cr10 resembles that of P. subiasi. When compared to all Spanish Cryobius, the external morphological characters of Cr9.2 and Cr10 are also most consistent with the ones described for P. subiasi. However, some do not match, such as the body length, the size of the eyes, characters of the pronotum and the elytra (Table 4). Additionally, the collection sites of these two specimens are not matching the currently known distribution of P. subiasi. Morphology suggests both limited relatedness to P. subiasi and to P. cantabricus.

The species from the Alps – P. apenninus (Dejean, 1831) (Apennine Alps), P. subsinuatus (Dejean, 1828) and P. unctulatus (Austrian Alps) – do not form a clade in the combined– or the CO1 phylogeny (Table 3: “Alpine Cryobius”).

The sequences for the Nearctic Cryobius species were obtained from Genbank, most are CO1 except for one sequence of 28S for one specimen (P. riparius). A clade with all those species together was never recovered, neither in the combined tree, nor in the CO1-only tree. According to the combined phylogeny, the Nearctic species are divided into two separate but well supported groups, respectively, which are intermixed with Palearctic species.

One group includes two Canadian species (P. brevicornis + P. empetricola) and one from Alaska (P. nivalis). Within that group, the two Canadian species form a well-supported clade (BV = 100, BP = 100). The other Nearctic species are grouped with the Japanese P. kurosawai Tanaka, 1958 which is closest to P. riparius (BV = 79, not recovered in the Bayesian inference). The relative position of these groups differs in the CO1 phylogeny, but the supports are lower.

Pterostichus (Cryobius) cf. anatolicus Jedlička, 1963 from Northeast Turkey is placed at the base of Cryobius as a sister to all other species, but without support. In addition to the combined phylogeny, this basal position of P. cf. anatolicus is also recovered by the analyses of CO1 and 28S only. However, in all cases this position is not well supported. 16S and 18S were not sequenced for this specimen.

This study is the first to provide phylogenetic data on the Pterostichus subgenus Cryobius. However, Cryobius was not tested for monophyly, but first molecular support for a monophyly of Cryobius is provided. Subsequent studies with a more comprehensive sampling of the genus Pterostichus are needed to further address this issue. Such subsequent studies would be important to verify the current phylogenetic status of Cryobius which is based on morphological clues (Ball 1966).

Within the phylogeny of Cryobius, no distinct clades of the former Haptoderus s.str. or Pyreneorites were recovered. Although not all species originally assigned to these two taxa were included in this study, the type species for Haptoderus (P. pumilio) and Pyreneorites (P. pusillus) were included.

In his revision of the genus Haptoderus, Jeannel (1937) distinguished the subgenera Haptoderus s.str. and Pyreneorites by the punctuation of the metepisterna, lacking in Pyreneorites, a character not mentioned by Ball (1966). Regarding the specimens of this study (Fig. 4), this character was reliable except for P. pumilio, for which the punctuation was sparse or absent in many specimens. Only the individual Cr27 showed a clear punctuation. Considering that this character is not mentioned in current determination keys, it could be that the punctuation of the metepisterna is a varying and hence unreliable character within Cryobius. Furthermore, Jeannel (1937) only used morphological clues for his separation of the former three subgenera. This led to the fact that the species of Pyreneorites were not the only ones exclusively distributed in the Pyrenees, as several species of Haptoderus s.str. were, according to Jeannel, also restricted to the Pyrenees (Jeannel 1937). In light of this, the synonymy of Pyreneorites with Haptoderus by Ball (1966) appears reasonable.

Ultimately, the synonymy of Haptoderus and Cryobius by Bousquet (1999) is also supported by the molecular data obtained in this study, as the specimens of Haptoderus (incl. Pyreneorites) and Cryobius form a monophyletic group.

The focus of this study was on the Pyrenean and Cantabrian Cryobius species. The results show that there are three lineages in the Pyreneo-Cantabrian mountain massifs, (i) one made up by P. (C.) pumilio-infimus, (ii) one by members of the P. (C.) cantabricus group and (iii) that formed by species of the P. (C.) abaxoides-amoenus group. The addition of molecular data from other taxa inhabiting either the Pyrenees (e.g. P. (C.) amaroides (Dejean, 1828), P. (C.) amblypterus (Chaudoir 1868)) or the Cantabrian Mountains (P. (C.) ehlersi (Heyden, 1881)) may even show the existence of new lineages. It is expected that a complete taxon sampling will show the existence of a large monophyletic clade comprising all taxa from the Pyrenees and the Cantabrian Mountains, including the widely distributed species P. (C.) pumilio, which likely became secondarily adapted to montane and lowland forests of central Europe. Within this large clade, others are expected to be found including (i) taxa restricted to the Pyrenees or the Cantabrian Mountains, or (ii) others occupying both mountain systems as is the case of the P. (C.) cantabricus clade. This hypothesis agrees with that formulated by Ortuño and Zaballos (1992) about a progressive colonization of the north Iberian mountains by ancestral Cryobius coming from the Pyrenees, perhaps from the onset of the Pleistocene.

In the combined tree, the Pyrenean species P. infimus is grouped with the widely distributed P. pumilio. This arrangement is not recovered in the CO1, 28S and 18S single gene analyses. At species level, P. infimus is morphologically well characterized. Therefore, a false determination of the specimen is unlikely. Pterostichus infimus might be the sister species of P. pumilio. This would have to be investigated further by adding more genetic data of other specimens and species. Pterostichus infimus comprises three subspecies with unreliable morphological characters (Pupier 2011; Bousquet 2017). Hence, a genetic analysis of those populations could help to clarify the status.

The other Pyrenean Cryobius form two sister clades. Within those, the clades of P. abaxoides and P. pusillus, each represented by two specimens, are well supported. For both species and most Pyrenean taxa, several subspecies have been described, as expected for alpine beetle populations with low dispersal power and reduced gene flow. The use of barcoding might help to assign specimens to described subspecies or to test the validity of these subspecies. Additionally, a broader sampling is necessary – including the type localities of all the described taxa – to precise the distribution of those populations. A question to answer would be if those subspecies form isolated populations across the Pyrenees or whether there is genetic exchange, in which case they are probably not valid and should be synonymized.

The Cantabrian species form a supported clade in which the position of the two P. aralarensis subspecies is unexpected as they are placed separately. The genetic differentiation between these two subspecies could be explained by the large geographic distance between the sampling locations. One corresponds to the occidental Pyrenees (P. (C.) a. aralarensis)), the other to the Cantabrian Massif (P. (C.) a. asturicus)), respectively. In this case, it seems that the morphological differentiation has occurred at a slower rate than the molecular one. This hypothesis deserves further in-depth investigation.

According to Pupier (2011), Jeanne reported P. aralarensis from the Pic d’Anie (Eastern Pyrenees) being far East from its Eastern distribution border in Spain (Aralar Range). However, this record is dubious (Coulon and Pupier 2014, Serrano pers. com.)

The two P. cf. subiasi specimens (Cr9.2 and Cr10), share several ambiguous clues compared to what is known about this species. Regarding the morphological characters, the larger size of the specimens, the form of the fore angles of the pronotum, and the lack of setae at the base of the second elytral stria in both specimens are the most noticeable. Also, the collection sites of the two specimens are about 200 km apart from the two only known distribution areas of P. subiasi in Northwest Spain: the Sierra de Los Ancares (Lugo) and the Sierra del lnvernadero (Orense) (Serrano 2013). However, this species was described relatively recently (1992) compared to other Cryobius species. It is therefore likely that the current information about the distribution of P. subiasi is incomplete due to a lack of sampling. Still, the morphological ambiguity does not allow a clear assignment of Cr9.2 and Cr10 to P. subiasi. It could be that these specimens belong to an undescribed species possibly close to P. subiasi. To clarify this, more sampling in the respective distribution areas and, most importantly, more genetic information would be needed.

Pterostichus pumilio comprises two subspecies: P. pumilio pumilio and P. pumilio nevadensis (Jeannel, 1947). However, P. pumilio nevadensis is not considered in this study, as it is only described from the Sierra Nevada in southern Spain, a record which was questioned by Serrano (2013).

A monophyletic clade for P. pumilio was only recovered in the CO1 tree, but with low support. This is due to the close relation with the P. infimus specimen that is discussed in paragraph 4.3. above. Still, two groups of P. pumilio are recovered in the phylogeny: “pumilio excl. Pyrenees/Cantabria” and the “Pyreneo-Cantabrian pumilio”. Interestingly, Jeannel (1949) discussed the presence of two forms for P. pumilio. First, the Central European form, which was described as spadiceus (Dejean, 1828). This form is distributed at lower altitudes in France but can be found above 1500 m in the Alps and the Massif Central. Second, the pumilio s.str. form which is typical for the Pyrenees (no mention of the Cantabrian Range) where it reaches very high altitudes. It is characterized by a minor size and differences in the shape of the pronotum and elytral striae. Given the collection localities, the two P. pumilio groups shown by the combined tree could reflect those two population groups. However, a morphological analysis of the nine P. pumilio specimens collected for this study could not confirm most of the characters described for pumilio and spadiceus. The absence of the parascutellar stria in the pumilio s.str. form could be confirmed for the Pyrenean specimens whereas the specimens from the Massif Central (Cr.20 and Cr27.2) showed slight remains of this stria, as described for the spadiceus form. However this finding must be viewed with caution since Español and Mateu (1945) generally considered the parascutellar stria a non-reliable character.

The altitude of the collection localities between the two groups did not differ significantly. More sampling and sequencing of those two groups is needed to further investigate Jeannel’s hypothesis of two P. pumilio forms. The fact that the Cantabrian specimen Cr24 of P. (C.) pumilio is closely related to Pyrenean ones (Fig. 3) confirms that this species is an exception within Cryobius of these mountains, due to its dispersal and colonization abilities denoted by its wide distribution area. It should be noted that Cr24 also lacked the parascutellar stria, as found in Pyrenean specimens.

The position of three species with an Alpine to Eastern European distribution (P. apenninus, P. subsinuatus and P. unctulatus) in the phylogeny does not allow statements concerning the lineage, as the supports are low and the sampling is scarce. Though, the lack of molecular data on species from the Alps and other European mountain systems is a major limitation to investigate whether the Cantabrian Range and the Pyrenees were colonized via the Alps or vice versa.

The Nearctic Cryobius species did not form a monophyletic group but were intermixed with Palearctic species. The reason for this could be that these specimens were only based on CO1, except for the P. riparius chimera combining a CO1– and a 28S sequence (Table 3). It can be expected that an inclusion of more taxa, and genes with different levels of conservation, would change the arrangement of the Nearctic Cryobius. Yet, most of the publicly available sequences for Nearctic Cryobius are currently from CO1.

An interesting result is the grouping of two Nearctic species with the Japanese P. kurosawai (Hokkaido). Unfortunately, this specimen was represented by a 28S fragment only. In the combined tree, as well as in the 28S tree, it is closest to the North American P. riparius, a species that, according to Bousquet (2012), “ranges from central Alaska to eastern Alberta” (Canada). According to Morita (2002), P. kurosawai is reported from Western Russia (Primorskij Territory, Sakhalin Island) and Northern Japan (Hokkaido, Rishiri Island). This relation between Palearctic and Nearctic Cryobius is in line with several studies showing a disjunct distribution of Eastern Asian and Northern American flora and fauna (Kruckeberg 1983; Ball and Currie 1997; Zhou et al. 2012; Weng et al. 2016; Liu et al. 2017; Haas et al. 2020; Sugawara et al. 2021). This biogeographical pattern is often explained by migration events across the Bering land bridge between the Cretaceous and the Quaternary (Sanmartín et al. 2001). In fact, in his revision of Cryobius, Ball (1966) stated that the subgenus Cryobius would have no close relatives in North America but was morphologically similar to several Palearctic subgenera of Pterostichus. Therefore, he assumed that the North American Cryobius would have a Palearctic origin. Ball and Currie (1997) listed Cryobius species according to their geographic distribution assigning several species as either Beringian, Palearctic-Beringian or Nearctic-Beringian. Here, in contrast to Bousquet (2012), P. riparius is listed as a Palearctic-Beringian species.