Between 2005 and 2012, different collecting campaigns were organised on mountainous reliefs of the Orobie Alps where the species was already known to be present, viz. Alben, northern Grigna, Presolana, and Arera (Burlini 1955). In addition, in order to investigate the distribution of C. barii, additional collecting campaigns were performed on mountainous reliefs surrounding the previous massifs and suitable for the presence of the species (i.e., elevation above 1800 meters and consisting of limestone) (Figure 1). Collected individuals were labelled with geographic coordinates, date and host plant on which they were found; then preserved in absolute ethanol and stored at –20 °C. With the exception of the Corna Grande population, where only two individuals were collected, seven to ten individuals from each population were selected for the DNA extraction. Total genomic DNA was extracted from each individual through the non-destructive procedure described in Montagna et al. (2013) and purified using the Qiagen DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) following manufacturer’s instructions. The extracted DNA was quantified using NanoDrop2000 spectrophotometer and used as a template for PCR reactions. A fragment of the mitochondrial cytochrome c oxidase subunit I (COI), of approximately 730 bp, was amplified with primer C1-J-2183 / TL2-N-3014 (Simon et al. 1994). PCR were performed in 25 µL reaction mix containing: 1´ GoTaq reaction Buffer (10 mM Tris-HCl at pH 8.3, 50 mM KCl and 1.5 mM MgCl₂), 0.2 mM of each deoxynucleotide triphosphate, 0.5 pmol of each primer, 0.6 U of GoTaq DNA Polymerase and 30 ng of DNA. PCR conditions used the following thermal cycle parameters: 3 min at 95 °C followed by 35 cycles of 30 s denaturation at 95 °C, 30 s annealing at 50 °C and 1 min 20 s extension at 72 °C, with a final single extra extension step of 10 min at 72 °C. PCR products were directly sequenced, in both strands, with ABI 3100 Automated Capillary DNA Sequencer (Applied Biosystem, Foster City, CA, USA). The obtained electropherograms were edited and assembled in a consensus sequence using Geneious Pro 5.3. The consensus sequences were deposited in GenBank with accession numbers MK492325-MK492374.

Figure 1. 

Cryptocephalus barii Burlini, 1948 distribution. A) Geographic location of the Orobie Alps and Cryptocephalus barii distribution. Yellow dots indicate localities where the species was observed, red dots indicate localities investigated with extensive sampling campaigns in which the species was absent (the source map was downloaded from http://www.geoportale.regione.lombardia.it/ and elaborated with QGIS 3.4.1). B) Cryptocephalus barii picture acquired using a Canon 450D camera; the multilayered micrographs were processed with Zerene Stacker (Richland, WA, USA).

The obtained 53 COI sequences were aligned using MUSCLE (Edgar 2004) with default parameters. The alignment was checked for reading frame errors and termination codons with MEGA 7 (Kumar et al. 2016). Intra-population and inter-population nucleotide p-distances were calculated with MEGA 7 (Kumar et al. 2016). In order to evaluate correlation between nucleotide and geographic distance matrices, the geographic distances between sampling locations were computed using R software (R Core Team 2016) starting from the latitude and longitude coordinates. Matrices correlation was estimated by Mantel test (Mantel 1967), as implemented in ade4 R package (Dray and Dufour 2007).

Phylogeographic relationships within C. barii were investigated based on Bayesian inference using the software BEAST 2.5.1. (Bouckaert et al. 2014) with the bModelTest module (Bouckaert and Drummond 2017) for the evaluation of the substitution model. In addition, the dataset was tested for strict clock model against the non-clock model using a Bayes factor comparison. The marginal likelihoods of the two models were estimated by the stepping stone method implemented in MrBayes 3.2 (Ronquist et al., 2012), and then compared. According to the criterion reported in Kass and Raftery (1995), the strict clock model was preferred since the difference between the marginal likelihoods of the two models was not significant (p-value > 0.2). We used the alignment partitioned by codon positions, as suggested by preliminary analysis performed with Partition Finder 2.1.1 (Lanfear et al. 2016), and the strict clock model with a substitution rate set to 0.0177 ± 0.0003 My^-1 (Papadopoulou et al. 2010), previously used for congeneric species (Montagna et al. 2016).

The tree prior was set using the Constant coalescent Kingman model (Kingman 1982). The analysis was carried out using a random starting tree, running two Markov chains for 100 million generations and sampling every 10,000 generations. Finally, the same analysis was performed sampling from priors only to evaluate the priors that we applied to the various parameters. Convergence was evaluated with Tracer version 1.7.1 (Rambaut et al. 2018), and the two chains were combined with LogCombiner (Drummond et al. 2012), discarding 20% of the trees as burn-in. The combined set of trees was summarised as a maximum clade credibility tree with TreeAnnotator (Drummond et al. 2012).

In order to confirm the rooting position of the tree, we carried out a preliminary phylogenetic reconstruction using an alignment consisting of the haplotype sequences of C. barii and orthologous sequences of five species mined from GenBank (i.e., Cryptocephalus cristula Dufour, 1843, Cryptocephalus asturiensis Heyden, 1870, Cryptocephalus flavipes Fabricius, 1781, Cryptocephalus azurescens Escalera, 1914 and Pachybrachis sp.; accession numbers: HE600320, HE600302, KJ765877, HE600310, HF947529) used as outgroups. In this phylogenetic reconstruction, performed on a dataset of 15 sequences, the tree prior was set using the Yule model (Yule 1925; Gernhard 2008), while other settings were the same of the previous analysis. In order to take into account for introgression and hybridisation phenomena a network of haplotypes was inferred through the minimum spanning network method as implemented in the software PopART 1.7 (Leigh and Bryant 2015).

We used the coalescent-based program MIGRATE-N 3.7 (Beerli and Palczewski 2010), which estimate past migration rates between populations, to test whether C. barii species survived ice ages in situ on ice-free nunataks or on large mountain of refuge at the periphery. It would be obviously impossible to evaluate every possible model of migration so we chose a small set of models in order to test our hypothesis. The following six migration models have been evaluated (Suppl. matreial 1: Figure S1): (A) directional migration from west to east, assuming a colonisation from the west; (B) directional migration from east to west, assuming a colonisation from the east; (C) directional migration from south to north, assuming a colonisation from southern refugia; (D) directional migration from north to south, assuming a colonisation from a possible northern refugium (Pegherolo); (E) a mixed model with colonisation from southern (Grigna) and northern (Presolana) refugia (Lohse et al. 2011); (F) a model, suggested by the phylogenetic analysis, with directional migration from Grigna to Alben, then to Pegherolo and Corna Grande and then to the others (Arera, Concarena and Presolana). We estimated the mutation scaled effective population size θ = xN_eµ (x = 1, for mitochondrial DNA), where N_e is the effective population size and µ is the mutation rate, as well as mutation scaled migration rates M = m/µ, where m is the immigration rate per generation. We used the marginal likelihood values approximated with thermodynamic integration to compare models with natural logarithm Bayes factor and to calculate the probability of each model (P(model_i)=mL(model_i)/∑_jⁿmL(model_j); Beerli 2012).

We used the sequence model of Felsenstein (1984) and random starting genealogy. A preliminary analysis was run with parameter values inferred by an F_st-based method to obtain θ and M estimates that were used as initial values of that parameters in subsequent analysis. Prior distributions for θ and M were uniformly distributed with boundaries 0–0.1 and 0–50000, respectively. We performed four independent runs for each analysis, each consisting of a burn-in period of 25 million Markov chain Monte Carlo (MCMC) steps, followed by 100 million steps. Samples were recorded every 5,000 steps, resulting in a total of 20,000 recorded parameter values for each replicate. We used a static heating scheme with four chains with temperatures 1.00, 1.50, 3.00 and 1,000,000, in order to improve the estimation of marginal likelihood (Beerli 2009).

Since the suitable habitat of the species is currently between 1,800 and 2,100 meters of altitude, we hypothesise that during glacials the amount of areas with a suitable habitat increase, thus allowing the formation of corridors connecting the massifs where the species is present. The correspondence between the phylogeography of the species and possible corridors of suitable habitat connecting the different mountainous reliefs was evaluated building maps of the Orobie Alps highlighting areas above a certain altitude by QGIS 3.4.1 software (QGIS Development Team 2009). The average vertical thermal gradient used in this study is of 0.54–0.58 °C every 100 meters, as estimated for Alpine regions by Rolland (2003).

A dataset of 35 presence localities was generated from GPS-precision field-recorded points for the target species C. barii. Nineteen bioclimatic variables were downloaded from the web repository Worldclim.org (Hijmans et al. 2005) at 30” spatial resolution and cut to the extent of the European Alps (sensu Biondi et al. 2013) through the ‘Extract by Mask’ tool in ArcMap 10.0 (ESRI, 2010). After this process, variables were tested for possible multicollinearity through the ‘Band Collection Statistics’ tool in ArcMap 10.0 (ESRI, 2010), a correlation matrix was calculated and variables’ pairs exceeding the Pearson’s r value of 0.85 were discarded (Elith et al. 2006; Iannella et al. 2018a; Iannella et al. 2018b).

For the modelling process, the ‘biomod2’ package (Thuiller et al. 2016) was used in R environment (R Core Team 2016). In particular, models for current and future climatic conditions were calculated through different sets of variables. Considering that many Global Climate Models (GCMs) are available for future climatic conditions, four GCMs were used in this paper, namely the BCC-CSM-1 (Wu et al. 2014), CCSM4 (Gent et al. 2011), IPSL (Marti et al. 2010) and the MIROC-CHEM (Watanabe et al. 2011). In particular, two scenarios of different radiative forcing were selected to observe the possible differences in conditions of medium and very high increase of radiative forcing, namely the 4.5 and the 8.5, for 2070.

Each species distribution model obtained from the different GCMs was processed through the MEDI algorithm (Iannella et al. 2017), a recent technique used to weight-average different predictions in one single model, thus avoiding predictions from one GCM and/or giving equal weight to models with low performances (see below).

In biomod2, Generalized Linear Models (GLMs, set with type = “quadratic”, interaction level = 3), Multiple Adaptive Regression Splines (MARS, set with type = “quadratic”, interaction level = 3), Generalized Boosting Model, also known as Boosted Regression Trees (BRT, set with number of trees = 5000, interaction depth = 3, cross-validation folds = 10) and Maxent (MaxEnt, set with maximum iterations = 5000) were selected as single modelling techniques (Phillips et al. 2006). Model was calibrated in current climatic conditions with the BIOMOD_Modelling function; all models’ performance were evaluated through the True Skill Statistics (TSS) (Allouche et al. 2006) and the Area Under the Curve (AUC) of the receiver operating characteristic curve (Phillips et al. 2006), with the initial 80% of the occurrence dataset used to calibrate the model and the remaining 20% for the validation. Then, the BIOMOD_EnsembleModelling function was used to obtain Ensemble Model for the target species, with the ‘wmean’ (weighted mean) algorithm used to merge each single model based on the respective performance scores. The BIOMOD_EnsembleForecasting function was further used to project the calibrated model to future climatic scenarios (Thuiller et al. 2016). A Minimum Convex Polygon (MCP) was generated on the basis of the presence data; all cartographic and spatial processes were managed in ArcMap 10.0 (ESRI, 2010).

The extensive sampling campaigns within the Orobie Alps and in neighbouring massifs suitable for the presence of C. barii (according with the proposed criteria, Materials and Methods), led to the collection of 60 individuals on eight massifs: the already known Mount Alben, Pizzo Arera, Presolana, and northern Grigna, with the addition of the newly discovered southern Grigna (hereafter reported as Grigna in association with the geographical neighbour northern Grigna), Corna Grande, Pegherolo, and Concarena (Figure 1). The performed extensive collecting campaigns make us to likely exclude the presence of C. barii on Mount Legnone, Mount Venturosa, Pizzo dei Tre Signori, Pizzo del Diavolo, and Mount Grona (Figure 1), mountainous reliefs presenting habitat suitable for the species and geographically close to previously known populations. All individuals were collected from mid-July to the end of August at an altitude ranging from 1,601 (only two samples, most likely transported by wind from higher altitudes) to 2,100 meters a.s.l., feeding or mating on Heliantemum nummularium, Hieracium spp. or Telekia speciosissima. The habitat of collection consists of grassland dominated by Sesleria coerulea and Carex sempervirens attributable to habitat code 6170 Alpine and subalpine calcareous grasslands (European Community Habitat’s Directive, 92/43/EEC).

DNA was extracted and COI amplified from 53 C. barii individuals (Table 1). Intra-population and inter-population mean pairwise nucleotide p-distance were 0.00047 (sd = 0.00064) [0–0.0018], and 0.011 (sd = 0.0083) [0.00049–0.024], respectively (Table 2). The population with the highest value of intrapopulation nucleotide p-distance was that of Presolana (0.0018, se = 0.0010), while the lowest values are associated with Alben, Corna Grande, and Grigna (0) (Table 2). Concerning the inter-population p-distance, the highest values were recovered comparing Grigna with other populations (mean p-distance = 0.023, sd = 0.0012), achieving the maximum value of 0.024 (se = 0.0057) when Grigna is compared with Alben, Corna Grande, and Concarena (Table 2). Some individuals collected in Presolana and Arera shared the same COI haplotypes, while all the other populations showed private haplotypes (Figure 2A). Positive correlation between geographic and nucleotide distances was detected by the Mantel test (r = 0.64, p-value = 0.001).

Figure 2. 

Maximum clade credibility tree and Minimum-spanning haplotype network. A Minimum-spanning haplotype network. Each colour represents a Cryptocephalus barii population. Circles represent the different haplotypes; their diameter is proportional to the haplotypes abundance. B Maximum clade credibility tree. Horizontal blue bars represent 95% HPD age confidence intervals for the nodes. Under the main lineage nodes is reported the divergence time in thousands of years before present. Above the nodes Bayesian posterior probabilities > 0.8 are reported, black asterisks indicate Bayesian posterior probabilities values < 0.80. Vertical coloured bars on the right of the tree indicate monophyletic clades, the colours identify the C. barii populations. Under the tree, in blue, the estimated surface temperature for the last 800,000 years (Hansen et al. 2013).

Based on the performed coalescence analysis, almost all the individuals from the same mountainous massif clustered together in monophyletic groups and are supported by high values of Bayesian posterior probability (BPP > 0.85; Figure 2B). However, in some cases, individuals from different massifs grouped together (Figure 2B); in details this phenomenon occurred for Arera and Concarena (two individuals; BPP = 0.92), Arera and Presolana (eight individuals; BPP < 0.5), Presolana and Arera (four individuals; BPP = 0.99) (Figure 2B).

Concerning the estimation of the divergence time among populations, the most ancient split, represented by the separation of Grigna lineage (tree rooted on outgroups) and all the remaining lineages occurred ~724,600 years before present (BP) (95% High posterior density (HPD) 1,140,700–389,900 years BP; BPP = 1) in correspondence with a period of warm climate, probably during the Pastonian or the Günz-Mindel interglacials (Figure 2B). The subsequent split, occurred about 244,700 years BP (95% HPD 417,900–107,500 years BP; BPP = 1) during the Mindel-Riss interglacial, determined the isolation of the Alben population and the next one, even if supported by values of posterior probability < 0.5, corresponded to the separation of Corna Grande-Pegherolo populations from individuals of Arera, Concarena, and Presolana. This last split is dated at about 149,100 years BP (95% HPD 263,900–68,900 years BP) corresponding to the Riss-Würm period, during the same interglacial the two populations of Corna Grande and Pegherolo separated one from the other ~87,800 years BP (95% HPD 197,700–24,500 years BP). The remaining populations (Arera, Concarena, and Presolana) diverged during the period from 162,700 to 83,700 years BP.

In order to understand how C. barii populations became isolated on different mountainous reliefs, six possible migration models were formalised and tested (Suppl. material 1: Figure S1). Bayes factors, calculated as double the difference of natural log marginal likelihoods (LBF), between two competing models strongly supported the model (E) (LBFs > 5, Prob_{model E} = 0.908, Figure 3A). In this model, Grigna and Presolana are considered source population, according with the hypothesis that C. barii survived ice ages both on ice-free refugia in the inner core of its range and on large mountain of refuge at its southern periphery, with unidirectional gene flow directed toward other populations which instead exchange migrants bidirectionally.

Figure 3. 

Most likely migration model and altitudinal habitat maps reporting suitable area for the presence of the species during cold periods. A Most likely migration model. Arrows indicate unidirectional flows (in black) and bidirectional flows (in gray) between populations. B–F Suitable altitudinal habitat maps. In light grey are reported the areas suitable for the presence of the species above a certain altitudinal threshold; the yellow ellipses are schematic drawing of Cryptocephalus barii populations; dendrograms showing the divergence events among populations, according to the tree in Fig. 1, are superimposed on the maps. At the bottom left of each map are reported: the minimum elevation at which climatic conditions are suitable for the presence of Cryptocephalus barii (inferred from present knowledge) corresponding to the altitudinal threshold used to draw the suitable areas, and the corresponding estimates of temperature variation in respect to the present. Abbreviations: A.T. = altitudinal threshold; CON = Concarena; PRE = Presolana; ARE = Arera; ALB = Alben; GRI = Grigna; NGR = northern Grigna; SGR = southern Grigna; CG = Corna Grande; PEG = Pegherolo.

The suitable altitudinal habitat maps, showing the increase of the areas appropriate for the C. barii survival and the available corridors connecting the present populations due to the decrease of the temperature, almost perfectly match with the topology achieved by the coalescent and the migration model analyses (Figure 3B–F). Interestingly, the populations inhabiting Arera, Presolana, and Concarena remained connected by habitat suitable for the species even when temperatures are only slightly lower than the present (Figure 3E). This last result is in agreement with the topology achieved by the phylogeographic analyses, where the relationships among these populations should be better represented by a polytomy or a network (BPP < 0.5; Figure 2).

Concerning the Species Distribution Modelling (SDM) analysis, multicollinearity among the nineteen bioclimatic predictors was prevented by discarding nine variables, keeping as predictors for the modelling process: the mean diurnal range (BIO2); the isothermality (BIO3); the maximum temperature of the warmest month (BIO5); the minimum temperature of the coldest month (BIO6); the mean temperature of the wettest quarter (BIO8); the mean temperature of the driest quarter (BIO9); the annual precipitation (BIO12); the precipitation of the wettest month (BIO13); the precipitation seasonality (BIO15); the precipitation of the driest quarter (BIO17). The corresponding correlation matrix is reported in Suppl. material 2: Table S1. Ensemble Models (EMs) resulted in very high scores of TSS (0.998) and AUC (0.988); the ‘wmean’ maps resulting show, for current climatic conditions, a predicted suitable area which narrowly encompasses the known distribution range (represented by a Minimum Convex Polygon, MCP) (Figure 4A). The suitable climatic conditions are, in fact, strictly limited within C. barii’s MCP, with few areas with low suitability in the surrounding mountains, and no suitable areas in the valleys (Figure 4A). The most contributive bioclimatic variables were resulted: BIO3 (34%), BIO17 (27%), and BIO2 (14%).

Figure 4. 

Predicted suitability for Cryptocephalus barii for current and future climatic conditions. Predicted suitability resulting from the Ensemble Modelling process performed over bioclimatic variables for Cryptocephalus barii, with the Minimum Convex Polygon built on the species’ presence sites for A current B 2070 – 4.5 scenario of radiative forcing, and C 2070 – 8.5 scenario of radiative forcing.

For future scenarios, the 2070_4.5 predictions resulted in a partial north-eastern shift of the habitat suitability, with an apparent increase of the compatible areas, which however show lower suitability with respect to the current situation (Figure 4B), even in the MCP area (Figure 5). For the 2070_8.5 scenario, a general decrease of suitability is observable, even though a range expansion is predicted; nevertheless, also in this scenario the suitability in the MCP is much lower than the current climatic Ensemble Model (Figure 4C, Figure 5).

Figure 5. 

Changes in habitat suitability for Cryptocephalus barii. Histogram reporting classes of habitat suitability calculated within the Minimum Convex Polygon built on Cryptocephalus barii presence sites through Ensemble Modelling process. Areas calculated for current and future climatic conditions (2070, 4.5 and 8.5 scenarios) are reported, respectively, in green, orange and red.