Baptiste Martinet
email: baptiste.martinet@umons.ac.be

Thomas Lecocq

Nicolas Brasero

Paolo Biella

Klára Urbanová

Irena Valterová

Maurizio Cornalba

Jan Ove Gjershaug

Denis Michez

Pierre Rasmont

Martinet Baptiste

Lecocq Thomas

Brasero Nicolas

Biella Paolo

Urbanová Klára

Valterová Irena
email: irena@uochb.cas.cz

Cornalba Maurizio

Ove

Michez Denis

Rasmont Pierre

Following the cold: geographical differentiation between interglacial refugia and speciation in the arcto-alpine species complex Bombus monticola (Hymenoptera: Apidae)

Keywords: Arcto-alpine regions, Bombus, cold-adapted species, integrative taxonomic approach, interrefugia differentiation, post-glacial warming

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Past climatic oscillations have led to significant changes in distributions of species.

However, species responses to climate change depend mainly on their eco-climatic requirements and tolerances [START_REF] Hewitt | The structure of biodiversity -Insights from molecular phylogeography[END_REF][START_REF] Stewart | Refugia revisited: Individualistic responses of species in space and time[END_REF]. Pleistocene and Quaternary climatic cycles triggered massive population movements resulting in periods of species range reductions (i.e. during cold periods when populations are restricted to refuge areas) for temperate species followed by periods of species range expansions (i.e. during warmer periods when populations recolonize at least portions of their initial range) [START_REF] Reinig | Die Holarktis. Ein Beitrag zur diluvialen und alluvialen Geschichte der Cirkumpolaren Faunen-und Florengebiete[END_REF]Hewitt, 2004a;[START_REF] Stewart | Refugia revisited: Individualistic responses of species in space and time[END_REF]. These population dynamics have fostered intraspecific divergence processes leading to differentiation and possibly speciation [START_REF] Avise | Phylogeography: The History and Formation of Species[END_REF][START_REF] Hewitt | The structure of biodiversity -Insights from molecular phylogeography[END_REF]. On the other hand, alternative demographic histories and subsequent differentiation patterns can be expected for cold-adapted species. Assessing accurately consequences of past climate change on differentiation and speciation process is a key element to better understand and predict the evolution of future biodiversity and to propose evidence-based mitigation strategies [START_REF] Rasmont | Climatic risk and distribution atlas of European bumblebees[END_REF].

While population dynamics of temperate species fostered by past climatic events and their consequences have been the focus of abundant research [START_REF] Zagwijn | The beginning of the Ice Age in Europe and its major subdivisions[END_REF][START_REF] Taberlet | Biodiversity at the intraspecific level: The comparative phylogeographic approach[END_REF][START_REF] Hewitt | Post-glacial re-colonization of European biota[END_REF][START_REF] Stewart | Refugia revisited: Individualistic responses of species in space and time[END_REF], cold-adapted species have received comparatively far less attention to date [START_REF] Mardulyn | Testing phylogeographic hypotheses in a euro-siberian cold-adapted leaf beetle with coalescent simulations[END_REF]. Contrary to temperate taxa, cold-adapted species are thought to have reached their largest distribution range during the Ice Ages [START_REF] Hewitt | Quaternary phylogeography: the roots of hybrid zones[END_REF].

The post-glacial warming and subsequent interglacial period is thought to have led to range contraction of such cold-adapted species toward the high latitude/altitude areas [START_REF] Barnes | Genetic structure and extinction of the woolly mammoth, Mammuthus primigenius[END_REF][START_REF] Fedorov | Comparative phylogeography and demographic history of the wood lemming (Myopus schisticolor): implications for late Quaternary history of the taiga species in Eurasia[END_REF][START_REF] Hewitt | Quaternary phylogeography: the roots of hybrid zones[END_REF]. Such a population dynamic scenario can explain current allopatric patterns of species distributed in the Arctic and in southern mountains (i.e. arcto-alpine species) acting as interglacial refugia. These taxa have their current distribution in the relicts (refugia) of a widespread distribution fragmented by post-glacial warming [START_REF] Reinig | Die Holarktis. Ein Beitrag zur diluvialen und alluvialen Geschichte der Cirkumpolaren Faunen-und Florengebiete[END_REF][START_REF] Mardulyn | Testing phylogeographic hypotheses in a euro-siberian cold-adapted leaf beetle with coalescent simulations[END_REF]Dellicour et al., 2014). In Europe, due to interglacial periods, arcto-alpine species exhibit a strong pattern of allopatry between Southern mountains (Pyrenees, Alps, Apennines, Balkans and Caucasus) and Northern areas (arctic regions of North Scandinavia and Russia). Such allopatric patterns have fostered and still foster gene flow disruptions leading to divergence and possibly speciation of cold-adapted species [START_REF] Avise | Phylogeography: The History and Formation of Species[END_REF][START_REF] Hewitt | Genetic consequences of climatic oscillations in the Quaternary[END_REF].

However, assessing species delimitation remains challenging because it requires to arbitrarily select variable traits whose accuracy continues to be debated [START_REF] Mayr | Systematics and the Origin of Species[END_REF][START_REF] De Queiroz | Species concepts and species delimitation[END_REF]Lecocq et al., 2015a). Moreover, it is quite difficult to comprehensively sample specimens for phylogeographic or speciation studies across vast inhospitable areas such as high altitude mountains and arctic areas [START_REF] Hewitt | Quaternary phylogeography: the roots of hybrid zones[END_REF]. This could lead to the underestimation of the variability within each allopatric population and to misunderstand the allopatric differentiation process.

The integrative taxonomy based on the unified species concept [START_REF] De Queiroz | Species concepts and species delimitation[END_REF] aims to overcome limitations due to unsettled adequacy of selected diagnostic traits and limited sampling. First, the approach considers multiple independent lines of evidence to evaluate inter-population differentiation processes and taxonomic statuses [START_REF] Schlick-Steiner | Integrative taxonomy: a multisource approach to exploring biodiversity[END_REF]Lecocq et al., 2015a,b). This reduces the likelihood of false taxonomic conclusions driven by single trait. Second, analyzing multiple traits to investigate inter-population differentiation allows to increase the amount of information available despite a limited sample size [START_REF] Lecocq | Molecular and chemical characters to evaluate species status of two cuckoo bumblebees: Bombus barbutellus and Bombus maxillosus (Hymenoptera, Apidae, Bombini)[END_REF].

Among potential organisms of interest for studying climatic oscillation consequences on cold-adapted species, bumblebees (Hymenoptera, Apidae, Bombus) represent a relevant biological system because some of them (i) live in the coldest areas inhabited by insects, (ii) have undergone diversification processes during the Pleistocene and Quaternary climatic cycles [START_REF] Hines | Historical biogeography, divergence times, and diversification patterns of bumble bees (Hymenoptera: Apidae: Bombus)[END_REF][START_REF] Duennes | Geographical patterns of genetic divergence in the widespread Mesoamerican bumble bee Bombus ephippiatus (Hymenoptera: Apidae)[END_REF]Martinet et al., 2015a;[START_REF] Rasmont | Climatic risk and distribution atlas of European bumblebees[END_REF][START_REF] Dellicour | Ecological niche modelling and coalescent simulations to explore the recent geographic range history of five widespread bumblebee species in Europe[END_REF]. Their interspecific and inter-population differentiations have been studied for a long time (e.g., Reinig, 1939). However, different diagnostic traits (morphological traits, DNA sequences, eco-chemical traits) have been used, resulting in conflicting biological conclusions (e.g. [START_REF] Gjershaug | Species status of Bombus monticola Smith (Hymenoptera: Apidae) supported by DNA barcoding[END_REF][START_REF] Williams | Genes suggest ancestral colour polymorphisms are shared across morphologically cryptic species in arctic bumblebees[END_REF]. Over the past few years, the efficiency of available diagnostic characters has been critically discussed and a merging of these traits in an integrative taxonomic framework has been proposed (e.g. Lecocq et al., 2015c). This provides the opportunity to efficiently delimitate species for a common cold-adapted bumblebee species with a strong pattern of allopatry.

Moreover, integrative taxonomy can help to define the subspecies status of allopatric populations (Lecocq et al., 2015a,c). In bumblebees, subspecies definition is traditionally based on color pattern variation but this diagnostic character requires an extensive overview of the inter-individual variability (Bertsch & Schweer, 2012a). However, color pattern has been shown to be unsuitable for taxonomic delimitation [START_REF] Vogt | Studien über das Artproblem. 1. Mitteilung. Über das Variieren der Hummeln. 1[END_REF]Bertsch & Schweer, 2012a;[START_REF] Carolan | Colour patterns do not diagnose species: quantitative evaluation of a DNA barcoded cryptic bumblebee complex[END_REF][START_REF] Williams | Genes suggest ancestral colour polymorphisms are shared across morphologically cryptic species in arctic bumblebees[END_REF] as well as for intraspecific variation study (Lecocq et al., 2015b,c).

Here, we investigated the potential inter-refugium differentiation and speciation within one of the most common arcto-alpine bumblebee species in European mountains (Rasmont et al., 2015): Bombus (Pyrobombus) monticola Smith, 1849. We sampled all the allopatric regions where the species is known (infra-specific taxa). We analyzed inter-population differentiation through multiple diagnostic traits: (i) a mtDNA marker (cytochrome oxidase I, COI), (ii) a nuDNA marker (phosphoenolpyruvate carboxykinase, PEPCK), and (iii) ecochemical traits (cephalic labial gland secretions, CLGS). Based on these traits, we developed an integrative taxonomic approach sensu Lecocq et al. (2015a,c) to assess the taxonomic status of major clades. In this approach, all taxonomic criteria used must be significantly differentiated to assign the species status.

MATERIAL & METHODS

Model Species

Bombus (Pyrobombus) monticola Smith, 1849 is an arcto-alpine species widespread in the alpine and sub-alpine stages of the highest mountain ranges of Europe with isolated populations in Northern Europe and Mediterranean mountains (Cantabrian Mountains, Pyrenees, Alps, Apennines, and Balkans but not Caucasus) [START_REF] Svensson | Pyrobombus lapponicus auct., in Europe recognized as two species: P. lapponicus (Fabricius 1793) and P. monticola (Smith, 1849) (Hymenoptera, Apoidea, Bombinae)[END_REF][START_REF] Kuhlmann | Checklist of the Western Palaearctic Bees (Hymenoptera: Apoidea: Anthophila)[END_REF][START_REF] Rasmont | Climatic risk and distribution atlas of European bumblebees[END_REF]. Bombus monticola was confirmed as an unique taxonomic unit by chemical (cephalic labial gland secretion and enzymology) and genetic analysis [START_REF] Svensson | Pyrobombus lapponicus auct., in Europe recognized as two species: P. lapponicus (Fabricius 1793) and P. monticola (Smith, 1849) (Hymenoptera, Apoidea, Bombinae)[END_REF][START_REF] Gjershaug | Species status of Bombus monticola Smith (Hymenoptera: Apidae) supported by DNA barcoding[END_REF] in comparison with its most similar taxon B. lapponicus (Fabricius, 1793). The analysis of [START_REF] Hines | Historical biogeography, divergence times, and diversification patterns of bumble bees (Hymenoptera: Apidae: Bombus)[END_REF], suggested that Bombus monticola diverged from its sister species B. lapponicus about 3 Ma. The species displays geographically differentiated color patterns [START_REF] Reinig | Die Verbreitungsgeschichte zweier fûr die Apenninen neuer boreoalpinen Hummelarten mit einem Versuch der Gliederung boreoalpiner Verbreitungsformen[END_REF] that have been used to define five phenotypically diagnosable allopatric subspecies (Tab.1 Fig 2, (Rasmont et al., 2015)): (i) B. monticola scandinavicus Friese, 1912 (Fennoscandia), (ii) B. monticola monticola Smith, 1849 (British Islands), (iii) B. monticola alpestris (= hypsophilus, Tkalcu, 1992) Vogt, 1909 (Alps, the Balkans, and the Olympus Mount), (iv) B. monticola rondoui Vogt, 1909 (Cantabrian Mountains and Pyrenees), and (v) B. monticola konradini Reinig, 1965 (Apennine Mountains) (Figs. 12). We define "monticola complex" as B. monticola ssp. + B. lapponicus and only "monticola" gathering exclusively all subspecies of B. monticola.

Sampling

We sampled 70 specimens including all B. monticola taxa (Appendix S1) from the entire known distribution area : B. monticola scandinavicus (n= 11) from North Scandinavia, B. monticola monticola (n=10) from the British Islands, B. monticola rondoui (n= 9) from the Pyrenees, B. monticola alpestris from the Alps (n=9), Balkans (n=3) and the Olympus Mountain (n=1) and B. monticola konradini (sensu [START_REF] Reinig | Die Verbreitungsgeschichte zweier fûr die Apenninen neuer boreoalpinen Hummelarten mit einem Versuch der Gliederung boreoalpiner Verbreitungsformen[END_REF] from the Central Apennines (Sibilini Mountains) (n= 5) and from the North Apennines (n=2). The North Apennines population, whose geographic distribution includes the highest peaks in the Apuan Alps, is separated by wide gaps not only from the Central Apennines populations but also from alpine alpestris (almost 230 km). We used the phylogenetically closely related species B.

(Pyrobombus) lapponicus (n=10) for comparison (see [START_REF] Cameron | A comprehensive phylogeny of the bumble bees (Bombus)[END_REF] and B. bimaculatus (Cresson, 1863) (n= 10) to root trees in our genetic analyses. All specimens were killed by freezing at -20°C. We considered all taxa without a priori taxonomic status and referred to them as scandinavicus, monticola, rondoui, alpestris, konradini, lapponicus, and bimaculatus (Tab.2). We further split konradini into konradini-N to indicate the Northern Apennines population and konradini-C to indicate the Central Apennines population.

Genetic Differentiation Analyses

To investigate the potential genetic differentiation between B. monticola taxa, we sequenced two genes that are commonly used in bee phylogenetic and phylogeographic studies (e.g. [START_REF] Pedersen | European bumblebees (Hymenoptera: Bombini) -Phylogenetic relationships inferred from DNA sequences[END_REF][START_REF] Cameron | A comprehensive phylogeny of the bumble bees (Bombus)[END_REF][START_REF] Williams | Unveiling cryptic species of the bumblebee subgenus Bombus s. str. worldwide with COI barcodes (Hymenoptera: Apidae)[END_REF], Dellicour et al., 2015): the mitochondrial gene COI and the nuclear gene PEPCK. We performed DNA extraction protocol, polymerase chain reaction (COI primers Apl2013/Aph2931, Pedersen 2002; PEPCK primers FHv4/RHv4, [START_REF] Cameron | A comprehensive phylogeny of the bumble bees (Bombus)[END_REF], sequencing procedures, and DNA sequence alignment using the method described in Lecocq et al. (2013). We uploaded the resulting COI (938 bp) and PEPCK (925 bp) sequences in GenBank (accession numbers Appendix S1).

We investigated the potential genetic differentiation within B. monticola through haplotype network analyses and phylogenetic inference. We carried out the analyses for each gene individually. We used the median-joining method to produce haplotype networks with Network 4. 6.1.0 (www.fluxus-engineering.com). We weighted transversions twice as high as transitions to reconstruct the network (Lecocq et al., 2015).

In phylogenetic analyses, we analyzed each gene with maximum parsimony (MP), maximum likelihood (ML) and Bayesian (MB) methods. We carried out maximum parsimony analyses (heuristic method) using Seaview 3.2 [START_REF] Galtier | SEAVIEW and PHYLO_WIN: two graphic tools for sequence alignment and molecular phylogeny[END_REF] with 1,000,000 replicas.

Only high quality trees and the majority rule 50% consensus tree were conserved. For ML and MB, each gene was partitioned as follows: (i) the nuclear gene (PEPCK) into two exons and two introns and (ii) each nuclear exon and (iii) the mitochondrial gene (COI) by base positions (first, second and third nucleotide) to define the best substitution model with JModelTest Server 2.0 [START_REF] Posada | jModelTest: Phylogenetic model averaging[END_REF] using the corrected Akaike information criterion. Bestfitting substitution models :( i) for COI : GTR+I (1st position), TIM2+I (2nd position), TrN+G (3rd position) ; (ii) for PEPCK intron 1: TPM1 uf +I ; (iii) for PEPCK exon 1: HKY+I (1st position), JC (2nd position), TrN+I (3rd position) ; (iv) for PEPCK intron 2: TrN+I ;(v) for PEPCK exon 2 : JC (1st position), JC (2nd position), JC (3rd position). For maximum likelihood analyses (ML), we performed 10 independent runs in GARLI 2.0 for both genes [START_REF] Zwickl | Genetic algorithm approaches for the phylogenetic analysis of large biological sequence datasets under the maximum likelihood criteria[END_REF]; the topology and-ln L was the same among replicates. Only the run with the highest likelihood was saved. We assessed statistical significance of nodes with 10,000 nonparametric bootstrap replicates. We considered a topology well supported (high confidence) whenever the bootstrap value (branch supports) was greater than 85% [START_REF] Hillis | An empirical test of bootstrapping as a method for assessing confidence in phylogenetic analysis[END_REF]. We carried out Bayesian inference analyses (MB) with MRBAYES 3.1.2 [START_REF] Ronquist | MrBayes 3: bayesian phylogenetic inference under mixed models[END_REF]. We achieved 10 independent analyses for each gene (100 million generations, four chains with mixed models, default priors, saving trees every 100 generations). Then we removed the first 10,000,000 generations as burn-in procedure. Then a majority-rule 50% consensus tree was constructed. Only branch supports (topologies) with high posterior probabilities (≥ 0.95) were considered statistically significant [START_REF] Wilcox | Phylogenetic relationships of the dwarf boas and a comparison of Bayesian and bootstrap measures of phylogenetic support[END_REF]. We (re) rooted all trees with the taxon B. bimaculatus.

To recognize species threshold, we used a Bayesian implementation of the general mixed Yule-coalescent model (bGMYC) for species delimitation based on COI tree (Reid & Carstens 2012, see an example of the use of the approach in Lecocq et al. 2015c). These analyses were performed with "bGMYC" R packages [START_REF] Reid | Phylogenetic estimation error can decrease the accuracy of species delimitation: a Bayesian implementation of the general mixed Yulecoalescent model[END_REF]. The stationarity and the modal coalescent/Yule ratio have been assessed to continue the analysis.

A range of probabilities > 0.95 was considered as strong evidence that taxa were conspecific while a range of probabilities < 0.05 suggested that taxa was heterospecific [START_REF] Reid | Phylogenetic estimation error can decrease the accuracy of species delimitation: a Bayesian implementation of the general mixed Yulecoalescent model[END_REF]. Since bGMYC required ultrametric trees, we performed a phylogenetic analysis with BEAST 1.7.4 (Drummond et al. 2012) using a phylogenetic clock model to generate a posterior distribution of trees (length of the MCMC chain: 1 billion generations). The first million sampled trees as burn-in, using the maximum clade credibility method and setting the posterior probability limit to 0. We based the bGMYC analysis on 1000 trees sampled every 10 000 generations. For each of these 1000 trees, the MCMC was made of 100 000 generations, discarding the first 90 000 as burn-in and sampling every 100 generations.

Posterior probability distribution have been applied against the first sample tree to provide a "heat' map".

Molecular Clock -Estimating Divergence Time

Following the approach of [START_REF] Duennes | Geographical patterns of genetic divergence in the widespread Mesoamerican bumble bee Bombus ephippiatus (Hymenoptera: Apidae)[END_REF] and Lecocq et al. (2013), we analyzed the COI dataset in BEAST v1.7.2 [START_REF] Drummond | BEAST: bayesian evolutionary analysis bysampling trees[END_REF] to estimate the divergence time among different clades. Using the GTR + I model selected by jModeltest, we ran Markov chain Monte Carlo simulations with the coalescent constant population size tree model and the relaxed clock model. Considering that no fossils of Pyrobombus species are available, the phylogeny is calibrated with a range date from a molecular study. We specified a range of possible substitution rates which includes the extreme rate for insect mitochondrial genes recorded in the literature (e.g. [START_REF] Duennes | Geographical patterns of genetic divergence in the widespread Mesoamerican bumble bee Bombus ephippiatus (Hymenoptera: Apidae)[END_REF] using a flat prior ranging from 1×10-9 to 1×10-7 substitutions site-1 and year-1. Simulations were run for 300 million generations, sampling every 1000 generations. Four independent runs were assessed in Tracer v1.4.1 (Rambaut & Drummond, 2016) to confirm convergence, determine burn-in, and examine the effective sample size of all posterior parameters. Log files from each run were combined in LogCombiner v1.6.1 (Rambaut & Drummond, 2016) for final parameter estimates.

Eco-Chemical Traits Differentiation

We focused on the most studied eco-chemical trait involved in the bumblebee premating recognition [START_REF] Baer | Bumblebees as model organisms to study male sexual selection in social insects[END_REF][START_REF] Ayasse | Chemical ecology of bumble bees[END_REF]: the cephalic labial gland secretions (CLGS). These secretions are complex mixtures of mainly aliphatic compounds synthesized de novo by male cephalic labial glands [START_REF] Coppée | Intraspecific variation of the cephalic labial gland secretions in Bombus terrestris (L.) (Hymenoptera: Apidae)[END_REF][START_REF] Lecocq | Molecular and chemical characters to evaluate species status of two cuckoo bumblebees: Bombus barbutellus and Bombus maxillosus (Hymenoptera, Apidae, Bombini)[END_REF]Žáček et al., 2013). We identified the main component as the compound that had the highest relative concentration (RA) among all compounds of CLGSs at least in one specimen of the taxon.

The CLGS are species-specific blends with some inter-population variations and are, subsequently, commonly used for species discrimination and assessment of intraspecific variability in bumblebees (review in Lecocq et al., 2015a). We extracted the CLGS with 400 µl of n-hexane, according to De [START_REF] De Meulemeester | A simplified protocol for bumble bee species identification by cephalic secretion analysis[END_REF] and [START_REF] Brasero | First chemical analysis and characterization of the male species-specific cephalic labial gland secretions of South American bumblebee[END_REF].

Samples were stored at -40° C prior to the analyses.

We qualified the CLGS composition of each sample by gas chromatography-mass spectrometry using a Focus GC (Thermo Scientific) with a non-polar DB-5 ms capillary column [5% phenyl (methyl) polysiloxane stationary phase; column length 30 m; inner diameter 0.25 mm; film thickness 0.25 μm] coupled to DSQ II quadrupol mass analyser (Thermo Scientific) with 70 eV electron impact ionization. We identified each compound using the retention times and mass spectra of each peak, in comparison to those from the National Institute of Standards and Technology library (NIST, U.S.A) database. We determined double bound positions (C=C) by dimethyl disulfide (DMDS) derivatization [START_REF] Vincenti | Determination of Double Bond Position in Di unsaturated Compounds by Mass Spectrometry of Dimethyl Disulfide Derivatives[END_REF].

We quantified the CLGS compounds with a gas chromatograph Shimadzu GC-2010 system (GC-FID) equipped with a non-polar SLB-5 ms capillary column [5% phenyl (methyl) polysiloxane stationary phase; column length 30 m; inner diameter 0.25 mm; film thickness 0.25 μm] and a flame ionization detector. We quantified the peak areas of compounds in GC solution postrun (Shimadzu Corporation) with automatic peak detection and noise measurement. The relative areas (RAs, expressed in %) of compounds in each sample were calculated by dividing the peak areas of compounds by the total area of all compounds. We excluded compounds for which RA were less than 0.1% for all specimens [START_REF] De Meulemeester | A simplified protocol for bumble bee species identification by cephalic secretion analysis[END_REF]. The data matrix for each taxon was based (Table S2) on the alignment of each relative proportion of compound between all samples performed with GCAligner 1.0 (Dellicour & Lecocq, 2013).

For GC/MS and GC-FID analyses, we injected 1 μL, using a splitless injection mode (injector temperature of 220° C) and helium as carrier gas (1 mL/min, constant velocity of 50cm/s). The oven temperature (of the column)was programmed isothermal at 70° C for 2 min and then rising from 70° C to 320° C at a rate of 10°C/min. The temperature was then held isothermal at 320° C for 5 min.

To facilitate the alignment of compounds and their identification, before each sample injection, a standard (Kovats) was injected containing a mix of hydrocarbons (alkanes) from C10 (decane) to C40 (tetracontane). Kovats indices were calculated with GCKovats 1.0 according to the method described by Dellicour & Lecocq (2013).

We performed statistical comparative analyses of the CLGSs using R environment (R Development Core Team, 2013) to detect CLGS differentiations between B. monticola taxa.

We used a clustering method, computed with the unweighted pair-group method with average linkage (UPGMA) based on Canberra distance matrices (RA of each compound) (R package ape; [START_REF] Legendre | Developments in environmental modelling 20[END_REF][START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF], to detect the divergence between taxa in the CLGS composition. We assessed the uncertainty in hierarchical cluster analysis using p-values calculated by multiscale bootstrap resampling with 100,000 bootstrap replications (significant branch supports >0.85) (R package pvclust, Suzuki & Shimodaira 2011). We assessed CLGS differentiations between taxa by performing a permutation multivariate analysis of variance using distance matrix (PerMANOVA) (R package vegan, [START_REF] Oksanen | Tertiary Vegan: Community Ecology Package[END_REF]. When a significant difference was detected, we performed a pairwise multiple comparison with an adjustment of p-values (Bonferroni correction) to avoid type I errors. We determined specific compounds of each taxon (indicator compounds) with the indicator-value (IndVal) method (Dufrêne & [START_REF] Dufrene | Species assemblages and indicator species: The need for a flexible asymmetrical approach[END_REF][START_REF] Claudet | Assessing the effects of Marine Protected Area (MPA) on a reef fish assemblage in a northwestern Mediterranean marine reserve: identifying community-based indicators[END_REF]. This value is the product of relative concentration and relative occurrence frequency of a compound within a group.

The statistical significance of an indicator compound (threshold of 0.01) was evaluated with a randomization procedure.

Morphological Analyses

To investigate diagnostic morphological character for species identification and new taxa description (not for species delimitation), a total of 60 workers were subjected to morphological analysis to discriminate B. lapponicus, B. konradini and B. monticola. We included only workers to have a sufficient sampling and because the differences in metabasitarsus measurements were more pronounced in females than males. We selected the maximum length and width metabasitarsus distance in reference to the work of [START_REF] Gjershaug | Species status of Bombus monticola Smith (Hymenoptera: Apidae) supported by DNA barcoding[END_REF] and we calculated the ratio (max length/max width) of these two measures to reduce the effect of body size on this morphological analysis. One picture was taken for each measurement and specimen using a binocular coupled with a digital camera (Nikon D70). The specimen was positioned in a way to maximize the focus of the metabasitarsus. The maximum metabasitarsus distance was measured on the picture with the software Imagej 1.5 [START_REF] References Abràmoff | Image Processing with Imagej[END_REF]) (Tab. 1, Appendix S3). Kruskal-Wallis analyses (Kruskal-Wallis test and multiple comparison test after Kruskal-Wallis, "pgirmess" R-package, [START_REF] Siegel | Non parametric statistics for the behavioural sciences[END_REF] were performed using R (R Development Core Team, 2013) to compare the different studied taxa.

Data Integration and Decision Framework

Assuming that species diagnosis and inter-population differentiation are more efficient in a multiple evidence based approach [START_REF] De Queiroz | Species concepts and species delimitation[END_REF][START_REF] Schlick-Steiner | Integrative taxonomy: a multisource approach to exploring biodiversity[END_REF], we proposed a species delimitation hypothesis according to our genetic and CLGS criteria based on the method performed by Lecocq et al.(2015a) derived from the approach established by [START_REF] Schlick-Steiner | Integrative taxonomy: a multisource approach to exploring biodiversity[END_REF] . In this method, all criteria used in the integrative approach must be convergent to assign the specific status. This strict approach could lead to underestimate the species differentiation but reduces the taxonomic inflation (Lecocq et al., 2015a) We assigned the species status to a taxon (with a high degree of certainty) if this taxon: (i) was genetically differentiated in all genetic markers (unique haplotype); (ii) constituted a monophyletic group with high branch support; (iii) was significantly differentiated in CLGS compositions (including IndVal indicator compounds, PerMANOVA, high bootstrap values > 0.85) (Lecocq et al., 2015a). We assigned the subspecies taxonomic status to phenotypically distinct allopatric populations with differentiation in some traits but not all in order to highlight those populations displaying such differentiation (originality) and to reduce the "underestimate's risk" of our strict approach to assign species status by naming them as a subspecies [START_REF] Hawlitschek | Island evolution and systematic revision of Comoran snakes: why and when subspecies still make sense[END_REF][START_REF] Ennen | Clinal variation or validation of a subspecies? A case study of the Graptemys nigrinoda complex (Testudines: Emydidae)[END_REF][START_REF] Lecocq | The alien's identity: consequences of taxonomic status for the international bumblebee trade regulations[END_REF].

Identification and Type Revision

The type series of Bombus lapponicus konradini [START_REF] Reinig | Die Verbreitungsgeschichte zweier fûr die Apenninen neuer boreoalpinen Hummelarten mit einem Versuch der Gliederung boreoalpiner Verbreitungsformen[END_REF] are presently at the Zoologische Staatssammlung München and have been revisited for this study. The identification of other studied taxa was checked with traditional identification keys such as [START_REF] Løken | Studies on Scandinavian bumble bees (Hymenoptera, Apidae)[END_REF] and [START_REF] Gjershaug | Species status of Bombus monticola Smith (Hymenoptera: Apidae) supported by DNA barcoding[END_REF].

RESULTS

Inter-Taxa Differentiation

Haplotype network analysis revealed six unique haplotypes for COI and two for PEPCK (Fig. 3) within B. monticola taxa complex. konradini-C was the only taxon displaying unique COI (6.8 % sequence difference from monticola alpestris and 5.3% from lapponicus) and PEPCK (0.7% sequence difference from monticola and 0.97% from lapponicus) haplotypes in the ingroup. All phylogenetic analyses (MP, ML and MB) of each single gene showed a similar topology with clades corresponding to haplotype groups found in the networks. Analyses showed strong support for all groups, but the position of konradini-C was variable in the clade in our phylogenetic analyses, and hence remains uncertain (Fig. 3).

Phylogenetic analyses on PEPCK showed two main lineages within "monticola" (Fig. 3): the central Apennines lineage (konradini-C, hereafter referred to simply as konradini) and the main lineage (all other taxa). COI-based trees resolved konradini as the sister group to the outgroup B. lapponicus rather than to other lineages of B. monticola. Among these last ones, COI phylogenetic trees underlined some geographic subgroups within "monticola" (Fig. 3): (i) the northern Apennine lineage of "monticola" (described hereafter as mathildis ssp. nov.), (ii) a western group including taxa from Pyrenees (rondoui) and Scotland (monticola); (iii) an eastern-northern group including specimens from Sweden (scandinavicus) and Alps + Balkans + Olympus Mt. (alpestris).

In comparison to the ML and MB analyses for COI data, the tree generated for bGMYC analysis displayed difference (not biologically significant) mainly in the branching of mathildis ssp. nov. As discussed in literature, these differences were probably due to the different parameters used in the BEAST 1.7.4 software to calculate the bGMYC model and because this pairwise matrix (heat map) was plotted against a sample tree ( Lecocq et al., 2015a;[START_REF] Barraclough | Diversification in sexual and asexual organisms[END_REF]. The bGMYC analysis (Fig. 4) highlighted several entities with low probabilities (<0.05) to be conspecific with the other ones. These results match with the same taxa recognized in the COI tree (MP, ML, MB analyses, Fig. 3). Overall, the bGMYC suggested the delimitation of four prospective species (p < 0.05) within the monticola complex (and the comparison group) as in Fig. 3 : (i) one group including all lapponicus (bGMYC conspecificity probabilities between individuals included in the group, p >0.98-1), (ii) a group with all konradini from Central Apennines (p >0.99-1), (iii) one group with all bimaculatus (p >0.98-1), (iv) all monticola subspecies (p > 0.13-0. The relative abundance of several compounds was different compared with the relative abundance in other taxa of B. monticola. The IndVal method highlighted several unique and diagnostic compounds of konradini (Tab. 3; i.e. ethyl tetradecenoate, ethyl tetradecanoate, hexadec-7-en-1-ol ethyl octadecadienoate, dotriacontane, ethyl octadec-9-enoate). In particular, konradini was characterized by ethyl octadec-9-enoate with a relative abundance of 8.28% while it had very low relative abundance in other subspecies (median 0.57%). The discrimination between konradini and other B. monticola taxa was supported by maximal bootstrap support values (100%) (Fig. 3). This differentiation was confirmed by statistical analysis (PerMANOVA F= 29.36 p-value < 0.05, between konradini and other B. monticola taxa).

Taxonomic Status

Species status was confirmed for the comparison group B. bimaculatus and B. lapponicus. According to the mtDNA and nuDNA divergence along with the CLGS composition differentiation (including main compounds) (Tab. 4), the species status was assigned to konradini (detailed information is given in Supporting information, Appendix S3).

Bombus konradini was originally described by [START_REF] Reinig | Die Verbreitungsgeschichte zweier fûr die Apenninen neuer boreoalpinen Hummelarten mit einem Versuch der Gliederung boreoalpiner Verbreitungsformen[END_REF] konradini nov. status) and the population from the Alps (B. monticola alpestris). Indeed, considering the slight differentiation in COI (0.53% of divergence from alpestris) and the strong divergence in coat color from alpestris (Tab.1, Fig. 3), the north Apennines population should have a new subspecies status: B. monticola mathildis (detailed information is given in Supporting information, Appendix S3).

Divergence Times among Clades

Based on the COI data, the divergence between B. konradini (Central Apennines) and B. lapponicus was estimated with a median of 0.79 Ma (min 0.25 -max 1.9 Ma) at the end of the Günz-Mindel inter-glacial period. The divergence time between the outgroup bimaculatus and the clade "monticola-lapponicus-konradini" was estimated with a median of 2.40 Ma (min 1.14 -max 3.88 Ma). In contrast, the divergence time between lapponicus-konradini and the clade "monticola" was estimated with a median of 2.30 Ma (min 1.23 -max 4.11 MA.

These last two divergence times correspond approximately to the onset of glaciation events and the formation of the Bering Strait. The other monticola subspecies have diverged recently with an estimated time of 40,000-18,000 (min 7,500 -max 548,000) years ago.

Morphological Analysis

Measurements of the ratio between the maximum length and width of the metabasitarsus show significant differences (Kruskal-Wallis multiple comparison chi²= 32.757; all p-values < 0.05) between lapponicus and monticola alpestris, monticola mathildis ssp. nov. The ratio is also significantly different between konradini and monticola ssp. but not between monticola alpestris and monticola mathildis (Fig. S2, Appendix S3). However, between konradini and lapponicus, although our results present a clear trend which highlight a larger ratio for konradini, there is no significant differentiation. According to these results, konradini appears as intermediate between monticola sensu stricto (large metabasitarsus ratio) and lapponicus (small metabasitarsus ratio). Diagnostic morphological characters are summarized in Tab. 1. Holotype (present designation): 1 male, labeled: 1) "Italy, Emilia Romagna, Reggio Emilia, Villa Minozzo, Mt Cusna, 44.283492 N 10.401028 E, 2057m, 05.VIII.2015, S/Scabiosa sp, Rec. M. Cornalba, BMAR0431"; 2) (on red paper) "Holotype"; 3) "det. B.

Impact of new taxa in zoological

Martinet 2016 Bombus (Pyrobombus) monticola mathildis Martinet, Cornalba & Rasmont" (Fig. 1).

Paratype: 2 males have been located, designated and labeled as paralectotypes: "Italy, Emilia Romagna, Reggio Emilia, Villa Minozzo, Mt Cusna, 44.283492 N 10.401028 E, 2057m, 05.VIII.2015, S/Scabiosa sp, Rec. M. Cornalba, BMAR0432" and "Italy, Emilia Romagna, Reggio Emilia, Villa Minozzo, Mt Cusna, 44.282288 N 10.401603 E, 2055m, 12.VIII.2015, S/Carduus carlinifolius, Rec. M. Cornalba, BMAR0433".

DISCUSSION

Inter-Population Differentiation of B. monticola

The concordance between genetic differentiation, geographic distribution, and CLGS divergence of populations suggests a strong intraspecific structure between the subspecies of monticola (Fig. 3). The western subspecies (B. monticola rondoui from Pyrenees and B. monticola monticola from British Islands), the north Apennines population (B. monticola mathildis ssp. nov.) and the eastern-northern subspecies (B. monticola scandinavicus from Sweden and B. monticola alpestris from Alps, Balkans and Olympus Mt) constitute, with the COI marker, five differentiated groups in three main lineages which diverged recently (about 40,000 -18,000 years ago based on molecular clock) during the Pleistocene/Quaternary (Fig. 3). This could explain the weak divergence of PEPCK marker between monticola subspecies (recent divergence) because nuclear genes have a lower mutation rate than mitochondrial genes [START_REF] Lunt | The insect Cytochrome Oxidase I gene: evolutionary patterns and conserved primers for phylogenetic studies[END_REF][START_REF] Trunz | Comprehensive phylogeny, biogeography and new classification of the diverse bee tribe Megachilini: Can we use DNA barcodes in phylogenies of large genera?[END_REF]. Such time-estimation matches with the start of the last post-glacial warming. So, it appears that the geographic pattern is most likely a consequence of allopatric differentiation and genetic drift triggered by a range fragmentation subsequent to the last post-glacial warming. We speculate that, at the beginning of the current interglacial time, taxa have found refuge in southern mountainous areas (Alps and southern peninsulas of Europe) and Northern Europe by contraction of their distribution areas or range shifting [START_REF] Hewitt | Post-glacial re-colonization of European biota[END_REF][START_REF] Hewitt | Inferring glacial refugia and historical immigrations with molecular phylogenies[END_REF][START_REF] Petit | Glacial refugia: hotspots but not melting pots of genetic diversity[END_REF][START_REF] Stewart | Refugia revisited: Individualistic responses of species in space and time[END_REF], similarly to the boreo-montane leaf beetle (Chrysomelidae), Gonioctena pallida [START_REF] Mardulyn | Testing phylogeographic hypotheses in a euro-siberian cold-adapted leaf beetle with coalescent simulations[END_REF]. The resulting allopatry has fostered mtDNA differentiation along with minor differentiation of chemical reproductive traits similarly to what has already been shown for insular populations of bumblebees (Lecocq et al., 2013(Lecocq et al., , 2015a)).

Despite their relative geographic isolation, all other B. monticola allopatric taxa previously recognized by [START_REF] Reinig | Die Verbreitungsgeschichte zweier fûr die Apenninen neuer boreoalpinen Hummelarten mit einem Versuch der Gliederung boreoalpiner Verbreitungsformen[END_REF] are considered as conspecific based on our diagnostic criteria with a low geographic genetic and phenotypic differentiation (decision framework Tab.4; Fig. 3). They shared overall the same CLGS composition (except some low relative concentration differences) and are characterized by only a slight genetic differentiation. These low differentiations, particularly in CGLS composition, can simply be explained by the short time of divergence due to geographical isolation and intraspecific variability [START_REF] Lecocq | Molecular and chemical characters to evaluate species status of two cuckoo bumblebees: Bombus barbutellus and Bombus maxillosus (Hymenoptera, Apidae, Bombini)[END_REF][START_REF] Lecocq | The alien's identity: consequences of taxonomic status for the international bumblebee trade regulations[END_REF]. Within Bombus monticola alpestris, the three sampled populations (Alps, Balkans and Olympus Mountain) are clearly consubspecific.

Mountain Top Speciation: Bombus konradini nov.status.

Contrary to the situation within B. monticola, B. konradini nov. status displays larger genetic and chemical traits differentiations (Fig. 3). Allopatry has most likely shaped the reproductive trait (CLGS) differentiation as observed in other species (Lecocq et al., 2013).

The strong genetic differentiation of B. konradini could be explained by an earlier divergence from the common ancestor with other B. monticola lineages, most likely temporally close to the B. monticola -B. lapponicus complex divergence. Indeed, based on genetic differences in the 16S gene, [START_REF] Hines | Historical biogeography, divergence times, and diversification patterns of bumble bees (Hymenoptera: Apidae: Bombus)[END_REF] suggested that B. lapponicus and B. monticola diverged from each other about 3 Ma. In temperate species, the post Ice-Age recolonization of territories by relict populations (from refugia), could lead to a new shuffling of the genetic pool by re-contacting of these populations without speciation [START_REF] Coyne | Speciation[END_REF][START_REF] Hewitt | The structure of biodiversity -Insights from molecular phylogeography[END_REF]. The modification of geographical range could trigger genetic and CLGS differentiation. Indeed, it has been shown that the reproductive traits including cephalic labial glands secretions can differentiate from both sides of the physical barriers that may exist between the refuge areas (Lecocq et al., 2013a). The case of the new species status of B. konradini lends strength to the hypothesis that for cold-adapted taxa, climatic oscillations (i.e. interglacial periods) have led to species differentiation in mountain refuges after geographical separation. Further phylogeographic and phylogenetic studies, based on larger sampling (including additional closely related species) and other genetic markers, are needed to accurately assess these hypotheses.

Our integrative taxonomic decision framework supported and confirmed the species status of B. monticola compared with its morphologically closely related species (B. lapponicus) [START_REF] Løken | Studies on Scandinavian bumble bees (Hymenoptera, Apidae)[END_REF][START_REF] Svensson | Pyrobombus lapponicus auct., in Europe recognized as two species: P. lapponicus (Fabricius 1793) and P. monticola (Smith, 1849) (Hymenoptera, Apoidea, Bombinae)[END_REF][START_REF] Cameron | A comprehensive phylogeny of the bumble bees (Bombus)[END_REF][START_REF] Gjershaug | Species status of Bombus monticola Smith (Hymenoptera: Apidae) supported by DNA barcoding[END_REF].

Our results also supported the species status of konradini which is endemic of high-altitudes (>1800 mt) of the Central Apennines (Manino et al., 2007) (Fig. 3, Tab. 4). Concerning the eco-chemical traits (CLGS), konradini differed from the other B. monticola taxa by lightweight compounds (volatile molecules) which could have a long distance attractive role [START_REF] Ayasse | Mating behavior and chemical communication in the order Hymenoptera[END_REF]. Therefore, the differentiation of these compounds may be a significant pre-mating reproductive barrier or may simply reflect divergence times and drift. Besides, according to the results of COI marker, konradini could be more closely related to B. lapponicus (Fig. 3) than B. monticola taxa as suggested in the original description of [START_REF] Reinig | Die Verbreitungsgeschichte zweier fûr die Apenninen neuer boreoalpinen Hummelarten mit einem Versuch der Gliederung boreoalpiner Verbreitungsformen[END_REF]. However, the phylogenetic position of konradini is not completely resolved because of the different tree topologies between COI and PEPCK results.

The species status of B. konradini suggests that inter-glacial periods can lead to species differentiation in mountain refugia in cold-adapted taxa. Unlike the populations of the Alps, Pyrenees and Balkans, where the interconnection and thus the possibility of exchanges and conspecificity are likely, the population of the central Apennines is much more isolated from other mountain chains with a possible endemic speciation (Martin-Bravo et al., 2010).

Several studies have shown the presence of endemic taxa in the Central Apennines (e.g. in amphibians Mattoccia et al., 2011, Canestrelli & Nascetti, 2008, 2012, in reptiles Nagy et al., 2002, Joger et al., 2007, in turtles Fritz et al., 2006, in plants Conti et al., 2005;[START_REF] Fuente | Sarcocornia hispanica (Chenopodiaceae), a new species from the Iberian Peninsula[END_REF][START_REF] Frattaroli | Seed germination capability of four endemic species in the Central Apennines (Italy): relationships with seed size[END_REF]in bumblebees Lecocq et al., 2013). For example, Lecocq et al., (2013a)

Conservation Remarks on B. monticola Complex and the Practice of Integrative Taxonomy

Considering all taxonomic criteria in our integrative approach (Fig. 3, Tab. 4), we propose to conserve the subspecies status for five monticola taxa [START_REF] Hawlitschek | Island evolution and systematic revision of Comoran snakes: why and when subspecies still make sense[END_REF]Lecocq et al., 2015a[START_REF] Dellicour | Ecological niche modelling and coalescent simulations to explore the recent geographic range history of five widespread bumblebee species in Europe[END_REF] North Apennines (formerly included by Reinig within konradini). Although the usefulness of subspecies status in bumblebees has been criticized and debated [START_REF] Ebach | How objective is a definition in the subspecies debate?[END_REF] during recent decades [START_REF] Williams | The bumble bees of the Kashmir Himalaya (Hymenoptera: Apidae, Bombini)[END_REF]Bertsch & Schweer, 2012a), we propose that these allopatric subspecies (partially isolated lineages) represent an important component and a useful pragmatic taxonomical unit for evolutionary biology and biological conservation in the evolutionary legacy of B. monticola (i.e. [START_REF] Waples | Evolutionary significant units and the conservation of biological diversity under the Endangered Species Act. Evolution and the Aquatic Ecosystem[END_REF][START_REF] Patten | Diagnosability versus mean differences of Sage Sparrow subspecies[END_REF][START_REF] Phillimore | Are subspecies useful in evolutionary and conservation biology[END_REF][START_REF] Rasmont | An overview of the Bombus terrestris (L. 1758) subspecies (Hymenoptera: Apidae)[END_REF][START_REF] Patten | Subspecies' and 'race' should not be used as synonyms[END_REF][START_REF] Crowhurst | Genetic relationships of hellbenders in the Ozark highlands of Missouri and conservation implications for the Ozark subspecies (Cryptobranchus alleganiensis bishopi)[END_REF][START_REF] Braby | The subspecies concept in butterflies: has its application in taxonomy and conservation biology outlived its usefulness[END_REF][START_REF] Sackett | Evidence for two subspecies of Gunnison's prairie dogs (Cynomys gunnisoni), and the general importance of the subspecies concept[END_REF]. These differentiations could be local adaptations to particular environments [START_REF] Avise | Phylogeography: The History and Formation of Species[END_REF][START_REF] Frankham | Introduction to Conservation Genetics[END_REF][START_REF] Braby | The subspecies concept in butterflies: has its application in taxonomy and conservation biology outlived its usefulness[END_REF]Lecocq et al., 2013a). Therefore, subspecies classification seems suitable to reflect the intraspecific differentiation within B. monticola taxa.

The monticola complex is a stunning example of the difficulty, in taxonomy, to define species or subspecies status of a population. Here the integrative taxonomy considering all criteria independently could provide strong pieces of evidence for decision. We assigned the subspecies taxonomic status to phenotypically distinct allopatric groups of populations with differentiation in some but not all criteria used in the integrative decision framework (i.e. conflict in selected criteria) [START_REF] Hawlitschek | Island evolution and systematic revision of Comoran snakes: why and when subspecies still make sense[END_REF][START_REF] Ennen | Clinal variation or validation of a subspecies? A case study of the Graptemys nigrinoda complex (Testudines: Emydidae)[END_REF]Lecocq et al., 2015a[START_REF] Dellicour | Ecological niche modelling and coalescent simulations to explore the recent geographic range history of five widespread bumblebee species in Europe[END_REF]. Taxonomical conclusions based only on the differentiation of one mitochondrial marker (e.g. COI barcoding) can lead to weak taxonomic hypotheses [START_REF] Andriollo | Two mitochondrial barcodes for one biological species: the case of European Kuhl's Pipistrelles (Chiroptera)[END_REF][START_REF] Mutanen | Species-Level Para-and Polyphyly in DNA Barcode Gene Trees: Strong Operational Bias in European Lepidoptera[END_REF]2016;[START_REF] Trunz | Comprehensive phylogeny, biogeography and new classification of the diverse bee tribe Megachilini: Can we use DNA barcodes in phylogenies of large genera?[END_REF] as mitochondrial differentiation may result from sex-specific characteristics, as lower dispersion for females [START_REF] Kraus | Male flight distance and population substructure in the bumblebee Bombus terrestris[END_REF][START_REF] Lepais | Estimation of bumblebee queen dispersal distances using sibship reconstruction method[END_REF], or mtDNA introgression or incomplete lineage sorting [START_REF] Bensasson | Mitochondrial Pseudogenes: evolution's miplaced witnesses[END_REF]Lecocq et al., 2015a). Taxonomic diagnosis based on multiple evidence (integrative taxonomy) is the best approach to avoid overestimation of species diversity which would lead to taxonomic inflation. Subspecies can be considered as a simple allopatric differentiation [START_REF] Mayr | Systematics and the Origin of Species[END_REF][START_REF] Patten | Null expectations in subspecies diagnosis[END_REF]. This procedure allows to assign a taxonomic status to any doubtful bumblebee taxa and points out these taxa for further taxonomic studies (Lecocq et al., 2015a). Moreover, despite the argument advanced by [START_REF] Williams | Genes suggest ancestral colour polymorphisms are shared across morphologically cryptic species in arctic bumblebees[END_REF], there is no case in bumblebees where the CLGS (mate recognition system) was not differentiated between two different species even when closely related bumblebee species have geographical distributions that do not overlap (e.g. B. terrestris (L.) and B. ignitus Smith, De Meulemeester et al., 2011 or between B. patagiatus and B. magnus where Bertsch & Schweer, (2012b) have shown the differentiation of their cephalic labial gland secretions).

Conservation Implication of the New Taxonomic Status of B. konradini nov. status

The new taxonomical status has implication for the red list assessments of the here studied European bumblebees according to the IUCN criteria [START_REF] Nieto | European red list of bees[END_REF]. While [START_REF] Rasmont | Climatic risk and distribution atlas of European bumblebees[END_REF] assess all taxa lumped into B. monticola, the new taxonomic status of B. konradini implies an evaluation of its conservation status independently from other B. monticola taxa. Bombus konradini was described as a rare, geographically very restricted taxon endemic to the central Apennines of Marche, Umbria, Lazio, Abruzzo and mostly occurring exclusively at elevations over 1800m [START_REF] Reinig | Die Verbreitungsgeschichte zweier fûr die Apenninen neuer boreoalpinen Hummelarten mit einem Versuch der Gliederung boreoalpiner Verbreitungsformen[END_REF]Ricciardelli & Piatti, 2003;[START_REF] Manino | Bumblebee (Bombus Latreille, 1802) distribution in high mountains and global warming[END_REF], Rasmont et al., 2015). The apparent scarcity of B. konradini could lead to significant genetic drifts (Ricciardelli & Piatti, 2003;[START_REF] Frankham | Introduction to Conservation Genetics[END_REF] that might significantly increase the species extinction risk [START_REF] Rasmont | Climatic risk and distribution atlas of European bumblebees[END_REF]. Indeed, according to [START_REF] Frankham | Introduction to Conservation Genetics[END_REF], small and isolated populations of a taxon are inherently more vulnerable to local extinction due to environmental and demographic stochasticity. It is therefore important to consider this new taxonomic status in our models and in our future backup plans (mitigation measures). monticola mathildis. With n= number of used specimens; * = significant differences (Kruskal-Wallis multiple comparison, p-value < 0.05). Tab.1. Range, conservation status and main morphological and colour pattern differences (male and female) between B. konradini nov. status, monticola subspecies including mathildis ssp. nov. and the similar species B. lapponicus according to [START_REF] Gjershaug | Species status of Bombus monticola Smith (Hymenoptera: Apidae) supported by DNA barcoding[END_REF][START_REF] Løken | Studies on Scandinavian bumble bees (Hymenoptera, Apidae)[END_REF]), Pittioni (1939) and original observations.

  95) including rondoui (IP >0.95-1), alpestris (p >0.95-1), scandinavicus (p >0.98-1), monticola (p >0.99-1) and mathildis ssp. nov. (p >0.99-1) which are significantly conspecific. The pairwise matrix (Fig. 4) shows more structure within B. monticola ssp. where the group displays more different haplotypes. These intermediate values of bGMYC (Fig. 4) between the different monticola lineages (genetic differentiation below the species differentiation threshold) are useful to discuss of subspecies concept. In chemical analyses, 103 compounds were detected; 82 in the CLGSs of B. monticola taxa (Appendix S2) except for konradini for which we detected only 50 compounds. The differentiation of CLGS composition between B. monticola taxa and outgroup species (B. lapponicus and B. bimaculatus) was conspicuous (IndVal; PerMANOVA F=115.63 and F= 122.52, p-value < 0.05; Fig.3). Except konradini, all other B. monticola taxa shared the same compounds with similar relative concentration (RA) (PerMANOVA F= 6.00 -13.20, p-value >0.05) (Appendix S2). Differences between konradini and other B. monticola taxa were particularly marked in the first half of the spectrum representing the most volatile molecules.

  as a subspecies typical of the northern and central Apennines, ranging from the provinces of Genova and Parma to L'Aquila. All other taxa were included in B. monticola but their color pattern (Tab.1) and/or differentiation in CLGS composition (minor quantitative differences) and/or in COI marker implied their assignation to the subspecies status. It is important to note the distinction of the north Apennines (province of Genoa and Parma to the provinces of Bologna and Lucca) monticola population (B. monticola mathildis ssp. nov) from the central Apennines taxon (B.

  nomenclature a) Bombus konradini sp. nov. status (more information in Appendix S3) Original taxonomic combination: Bombus lapponicus konradini Reinig 1965: 105 Locus typicus: Monti Sibillini, Central Apennine Mountains (Italy) Syntypes: 13 queens, 93 workers, 28 males Lectotype (present designation): 1 queen, labeled: 1)"Italia, Monti Sibillini, Nh. M. Vettore, Baumgrenze, 15-1600m, 14.6.61, Reinig".; 2) (on red paper) "LECTOTYPE; 3) "det. P. Rasmont 2015 Bombus (Pyrobombus) monticola konradini Reinig" (Fig.1). Paralectotype: 2 queens, 41 workers, 16 males have been located, designated and labeled as paralectotypes (Tab.2). In this series, only 1 queen (lectotype) and 21 workers from Marche, Umbria, Lazio and Abruzzo have been identified as Bombus konradini. The remaining paralectotypes (2 queens, 20 workers, 16 males) from Liguria, Emilia-Romagna and Toscana have been assumed as Bombus monticola mathildis. b) Bombus monticola mathildis Martinet, Cornalba & Rasmont ssp. nov. (more information in Appendix S3) Locus typicus: North Apennines, Emilia Romagna, Reggio Emilia, Villa Minozzo, Mt Cusna (Italy).

  provided evidence that the population of B. lapidarius (a Palearctic polytypic species) from the Southern Italian refugia has experienced genetic and CLGS differentiation during Quaternary glaciations leading to an incipient speciation process. Populations inhabiting the Mediterranean mountains, (e.g. the Apennines, one of the few mountain ranges in Europe arranged on a north-south axis) are characterized by a high genetic diversity (hotspot) with endemic taxa(Ruiz-La-Bourdette et al., 2012).The sympatry of two different species, defined by divergent taxonomic traits, reinforces the "species" status because individuals co-inhabit in the same area without hybridization. Our results suggest that B. monticola (sensu stricto) is absent in Central Apennines unless this could be due to a sampling bias. Such absence could result in lack of sympatry between B. monticola and B. konradini. Several hypotheses could explain the potential absence of B. monticola (sensu stricto) in Central Apennines: (i) for eco-climatic constraints, historical or competition reasons, this taxon has never inhabited this region or has disappeared, (ii) despite the significant observed differences (genetic, morphological and chemical traits), a limited hybridization between monticola and konradini could be still possible. Following this second hypothesis, along their contact zone, the subspecies mathildis, could represent an intermediate population resulting from some introgressions of the population living in Alps (alpestris). Our COI results suggest that the subspecies mathildis (low branch support) is closer to B. konradini than all other subspecies of monticola (sensu stricto) (Fig.3). However, phenotypic and chemical traits results do not tend towards this hypothesis (Figs.3C, S2). Although distinct, B. konradini could be the "replacement species" to B. monticola (sensu stricto) with similar eco-climatic constraints and filling the ecological niche in Apennines or a relict population of a near relative of B. lapponicus in Italy considering COI results. Additional ethological experiments (hybridization tests) and further genetic analyses (e.g. Microsatellite, SNPs, RAD-seq) are necessary to test these hypotheses of intermediate populations or replacement species in the context of taxonomic implications.

Fig. 1 .

 1 Fig. 1. Distribution map (Gall projection) of Bombus monticola (Rasmont & Iserbyt, 2014) and its traditional subspecies in Europe according to Rasmont (1983a). A, Bombus monticola scandinavicus queen, red area on the map; B, Bombus monticola monticola queen, dark green area; C, Bombus monticola rondoui queen, purple area; D, Bombus monticola alpestris queen, blue area; E, Bombus monticola mathildis ssp. nov Holotype male, pink area; F, Bombus konradini nov status Lectotype queen, yellow area. Light green dots indicate the occurrence of the taxon in the region. All photographs are by P. Rasmont.

Fig. 2 .

 2 Fig. 2. Photos of the studied bumblebee taxa. A, Bombus monticola scandinavicus queen; B, Bombus monticola monticola queen; C, Bombus monticola rondoui queen; D, Bombus monticola alpestris queen; E, Bombus monticola mathildis ssp. nov Holotype male; F, Bombus konradini nov status Lectotype queen. All photographs are by P. Rasmont.

Fig. 3 .

 3 Fig. 3. Genetic and chemical analyses within the monticola complex. (A) Majority rule (50%) consensus tree based on maximum likelihood analyses of COI. Values above tree branches are Parsimony bootstrap values/ maximum likelihood bootstrap values/Bayesian posterior probabilities. Only maximum likelihood and parsimony bootstrap values >70% and posterior probabilities >0.95 are shown. (B) Majority rule (50%) consensus tree based on maximum likelihood analyses of PEPCK. Values above tree branches are Parsimony bootstrap values/ maximum likelihood bootstrap values/Bayesian posterior probabilities. Only maximum likelihood and parsimony bootstrap values >70% and posterior probabilities >0.95 are shown. (C) 1) Dendrogram of cephalic labial gland secretion differentiation within monticola complex and B. bimaculatus. This cluster was obtained by hierarchical clustering using an unweighted pair-group method with arithmetic mean (UPGMA) based on a Canberra matrix calculated from the cephalic labial gland secretion matrix of B. bimaculatus (red), B. lapponicus (dark blue), B. konradini nov. status (green), B. m. rondoui (pink), B. m. scandinavicus (yellow), B. m. monticola (light blue), B. m. alpestris (black), B. m. mathildis ssp. nov. (orange). The values near nodes represent multiscale bootstrap resampling values (only values >80 of main groups are shown except nodes between B. monticola subspecies). 2) Principal component analysis (PCA) of cephalic labial gland secretion differentiation within monticola complex and B. bimaculatus: B. bimaculatus (red circles), B. lapponicus (dark blue circles), B. konradini nov. status (green circles), B. m. rondoui (pink circles), B. m. scandinavicus (yellow circles), B. m. monticola (light blue circles), B. m. alpestris (black

Fig 4 .

 4 Fig 4. Species recognition pairwise matrix. Species recognition pairwise matrix based on ultrametric tree of cytochrome oxidase 1 (COI) sequences with bGMYC pairwise probability of conspecificity plotted on a sample tree from BEAST. The colored matrix corresponds to the pairwise probabilities of conspecificity returned by the bGMYC method (color scale on the right of the figure). Black spots show the coalescent node for each species.

  the punctures on vertex is shiny, and there is a slight depression with some punctures near the compound the maximal width of the basitarsus is high (sensu[START_REF] Gjershaug | Species status of Bombus monticola Smith (Hymenoptera: Apidae) supported by DNA barcoding[END_REF]. The length of the metabasitarsus of these taxa is large (

  : B. monticola rondoui from Pyrenees, B. monticola monticola from British Islands, B. monticola scandinavicus from Fennoscandia and B. monticola alpestris from Alps, Balkans, Olympus Mountains and B. monticola mathildis ssp. nov. from

ACKNOWLEDGEMENTS

The authors thank the Abisko and Tarfala scientific stations (Sweden) for their welcome and their help in material collection. We acknowledge all people that helped us in our journey to the Abisko and Tarfala stations: P. Lakare and G. N. Rosqvist (University of Stockholm), M. Augner and L. Wanhatalo (Abisko Station), H. Savela (Oulu University, INTERACT project administration), J. Strand and T. Wikstrom (Lansstyrelsen i Norrbottens lan Naturvardsenheten, Lulea). The authors also thank the Parco Nazionale dei Monti Sibillini and the Parco Nazionale dell'Appennino Tosco-Emiliano for granting permission to collect in their respective territories. Special thanks go to P. Salvi (Sibillini), W. Reggioni (Appennino Tosco-Emiliano) and A. Cetkovic (University of Belgrade) for their help in the sampling. Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Belgian F.R.S.-FNRS. We thank also the two anonymous reviewers for their help to improve this manuscript. BM contributes as PhD student granted by Research Council of University of Mons and by the FRS-FNRS (Fonds de la Recherche Scientifique). PB contributes as PhD student funded by Czech Science Foundation (GAČR GP14-10035P) and by University of South Bohemia (GA JU 152/2016/P). Part of this work (I.V.) was supported by the Institute of Organic Chemistry and Biochemistry of the Academy of Sciences of the Czech Republic (project No. 61388963).

FUNDING

The research has received funding from the European Community's Seventh Framework Program, STEP Project (Status and Trends of European Pollinators, www.step-project.net, grant agreement no 244090, FP7/2007-2013) and from the INTERACT project.

AUTHOR CONTRIBUTIONS

Conceived