Most of the specimens used in this study were sampled during a specific research expedition carried out in July 2018 by two of the authors (M.M. and M.K.), with the addition of several previously collected specimens, mainly conserved in coll. Monguzzi. Specimens were collected usually at night, by searching on the limestone rocks near plants of Saxifraga caesia. Some more specimens were sifted from the same plants, and larvae were found in the roots of the saxifrage. Fresh imagoes and larvae were stored in vials with alcohol 95° immediately after collection and conserved at minus 23°C until they were processed. Specimens from previous collections had been conserved dry on card in entomological drawers and an attempt at DNA extraction from these specimens, not always successful, was made after a short rehydration.

Body length was measured in profile from anterior border of eyes to apex of elytra, excluding rostrum. Length/width ratios were measured from digital photographs, and were always taken at the maximum length and width of the respective parts in dorsal view. Genitalia were cleared with 10% KOH and carefully dissected; female genitalia and male genital sclerite were embedded in resin, male genitalia were mounted dry on the same card as the respective specimen. Photographs were taken using a Nikon P 6000 digital camera mounted on a Leica 6SE stereomicroscope, combining image stacks with Zerene Stacker. All images were cleaned and enhanced as necessary in Adobe Photoshop CS3.

Total DNA was extracted by placing the entire animal body in 400 μl of 5M guanidine-isothiocyanate, after separating the head + prothorax from the rest of the body, or cephalic capsule from the rest of the body for the larvae, to maximize DNA extraction. DNA extraction was destructive for the larvae; adult specimens were conserved as much as possible. Two different regions of the mitochondrial COI gene were amplified and independently analysed. One region was amplified with primers based on Hughes and Vogler (2004): fw C1-J-2183 (Jerry), 5’-CAACATTTATTTTGATTTTTTGG-3’ and rev L2-N-3014 (Pat), 5’-TCCAATGCACTAATCTGCCATATTA-3’. The other region was amplified with primers based on Folmer et al. (1994) modified as in Astrin and Stüben (2008): fw: LCO1490-JJ, 5’-CHACWAAYCATAAAGATATYGG-3’; rev: HCO2198-JJ, 5’-AWACTTCVGGRTGVCCAAARAATCA-3’. In addition, sequences corresponding to the Internal Transcribed Spacer 2 (ITS2) of the ribosomal DNA were amplified using primers based on Vahtera and Muona (2006): fw: 5’-GGGTCGATGAAGAACGCAGC-3’; rev: 5’-ATATGCTTAAATTCAGCGGG-3’. Amplification of DNA was done as follows: 15 min of initial denaturation (95°C) followed by 10 cycles of 30 sec at 94°C, 45 sec at 60°C to 50°C (lowering the annealing temperature in each cycle 1°C), 2 min at 72°C followed by 30 cycles of 30 sec at 94°C, 45 sec at 50°C, 2 min at 72°C and a final extension cycle of 15 min at 72°C. The reaction products were purified by agarose gel electrophoresis and successive purification from the gel. Sequencing was performed by an external service (Genechron, Roma). Both strands were sequenced. Forward and reverse chromatograms were checked with Chromas (http://technelysium.com.au/wp/chromas) using default parameters, and ambiguities were corrected manually. The COI sequences had no indels after alignment, no stop codons were detected and the translation to the amino acids was congruent with the protein sequence known for the other species of the genus. This was considered as a sufficient evidence to exclude NUMT pseudogenes. Multiple sequence alignment of both strands for both COI regions and ITS2 was performed with MEGA-X (Kumar et al. 2018), after reversing and complementing the reverse strand, with the Muscle alignment option. The COI sequences of the fragment amplified with the Jerry/Pat primers were aligned with those of other species of the genus used for inferring phylogeny of Dichotrachelus by Meregalli et al. (2018). The sequences of the fragment amplified with the Folmer primers were aligned with the few other Dichotrachelus Folmer sequences available, with the inclusion of one species of Cyclominae and one of Entiminae as outgroups. The ITS2 sequences were aligned with those of a few other species of Dichotrachelus kindly provided by C. Germann, with a species of the Molytinae genus Aclees as the outgroup. All the sequences were trimmed at the extremes to exclude the part corresponding to the primers and the final sequences were, respectively, 826 and 658 bp long for the two regions of COI and 702 bp long for ITS2 (including gaps for this last one, ITS2 634 bp long for the D. grignensis complex only). In the two regions of COI amplified, the first codon of the translated amino acid chain corresponded to the second site. In order to apply a model that allows for different substitution probabilities according to the position of the nucleotide in the codon, the first nucleotide was also trimmed and the analyses were therefore conducted on fragments of, respectively, 825 and 657 bp. All the sequences used for the phylogenetic analysis were deposited in GenBank.

Pairwise distance was calculated for the COI sequences with MEGA-X, implementing Tamura-Nei substitution model (Tamura and Nei 1993) with Γ-distributed rates between sites.

Bayesian Inference was estimated using MrBayes 3.2 (Ronquist et al. 2009). Two runs with 4 chains were run for 2 million generations, sampling every 500 generations. The chains were left free to sample all the models of the GTR family using reversible jump Monte Carlo Markov Chain (MCMC) (Huelsenbeck et al., 2004). Heterogeneity of substitution rates among different sites was modelled with a 4-categories discretized Γ distribution, with a proportion of invariable sites. For the two regions of COI, the matrix was partitioned so that substitution rates could vary according to the nucleotide position in the codon. A basic approach, with no constraints to the substitution rate for each site, was applied to ITS2. The first 25% generations were discarded (burn-in) and convergence was evaluated with the average standard deviation of split frequencies. Goodness of mixing was assessed looking at the acceptance rate of swaps between adjacent chains, following Ronquist et al. (2009). After a first analysis, the parameter “temperature” was set to 0.05 in order to improve swaps between chains for the COI sequences. The resulting consensus tree was examined with Figtree (Rambaut 2014).

An attempt to estimate the epoch of the splitting event between the two D. grignensis clades that resulted from the Bayesian inference, and within the Bergamasque Prealps populations, was carried out by applying a molecular clock to the results of the COI analyses. Criticism of molecular clocks has often been expressed, in particular because variation among taxa and intrinsic theoretical limits of the procedure have been evidenced (Bromham and Penny 2003; Schwartz and Maresca 2006; Gibbons 2012; etc.). However, the theoretical bases of the molecular clock have been thoroughly discussed (Ho and Lo 2013) and nowadays the Bayesian clock dating methodology has become a standard tool for estimating the timeline of the Tree of Life (dos Reis et al. 2015). Some peculiarities of Dichotrachelus must be considered. Firstly, it is impossible to calibrate the clock based on known age of independent events, let apart fossils, never recorded for these weevils, which would allow a more precise tree node dating. Moreover, the weevils’ generation time should be considered, since the overall rate of evolution is ultimately constrained by the turnover rate of individuals in populations, as reflected in generation time: species with lower generation turnover obviously tend to have lower substitution rates per site per million years (‘My’) (Thomas et al. 2010). Also the effect of metabolism may have influence (Gillooly et al. 2004), even though this last assumption has been criticized (Lanfear et al. 2007). The contrasting opinions underline uncertainties in the application of the molecular clock, which should be always taken with a margin of doubt. These Dichotrachelus have only a few months a year of active life, due to the harsh climate in their high altitude alpine habitats, and their life cycle seems to be at least two years long (Meregalli 1980). Various studies suggested that in Coleoptera the nucleotide substitution rate ranges around 0.015–0.018 substitutions per site per My (as, respectively, in Pons et al. 2010 and Papadopoulou et al. 2010), and this value was used here, but it cannot be excluded that indeed the true substitution rate is lower than 0.015 substitutions per site per My because of the metabolic constraints. As a consequence, this would increase the time of divergence between the phylogenetic lines.

Three different molecular clock models were tested with Bayes factor (B10): no-clock model, strict clock model, and uncorrelated gamma rate relaxed clock model (Lepage et al. 2007). MrBayes 3.2 was used to compute the marginal likelihood of the three models with the stepping stone algorithm (Xie et al. 2012): a 50 steps 2.5 million generations analysis sampling every 500 generations was performed using 2 runs with 4 chains, each under the GTR model. The Ln marginal likelihood for the no-clock, strict clock, and uncorrelated gamma rate model was respectively: = –2760, –2785, –2752 (825 bp fragment); –2642, –2622, –2612 (657 bp fragment). According to Kass and Raftery (1995) this is a very strong evidence in favour of the uncorrelated gamma rate relaxed clock, that we eventually used for dating.

Node ages were calculated with a rate of 0.018 nuclear substitutions per site per My per lineage, according to Papadopoulou et al. (2010). Uncertainty of the substitution rate is evaluated in Bayesian inference by changing the rate with a proposal mechanism at every generation. A Bayesian analysis with the same parameters of the phylogenetic analysis was performed (except the model and the parameters related to the clock).

The two taxa are extremely similar, and each population shows some morphological peculiarities; moreover, the very limited number of specimens available does not allow evaluation of the stability of the few distinctive characters; due to the paucity of specimens, no statistical analysis of possible biometric variation can be implemented. The following scheme was drawn from the specimens of D. orobicus from the type locality, Pizzo Camino, and surroundings, vs 8 specimens from the Grigna (Figs 2, 3).

The Bayesian inference gave phylogenetic trees with comparable topology for all the primer sets. All populations belonging to the D. grignensis complex clustered in a fully supported clade. Within this clade all the specimens from the Bergamasque Prealps formed a monophyletic group with respect to the specimens from the Grigna Mts with full, or almost full, support for both the regions of COI and the ITS2 (Figs 4, 5, 6).

According to our analysis, the two clades of the D. grignensis complex separated around 5 million years ago (‘Mya’), with similar results for both the 825 and 657 fragments; the median in-clade divergence among the populations of D. orobicus indicates that they became reciprocally isolated at most around 0.5 Mya, but probably in more recent times (Fig. 7).

The p-distance between the specimens of D. orobicus and those of D. grignensis is, for the 825 fragment, about 0.06 (6%), whereas the interpopulation p-distance for D. orobicus is always below 1%, varying between 0.004 and 0.008; for the 657 fragment the p-distance between D. orobicus and D. grignensis is 0.07 and among the various populations of D. orobicus it is lower than 0.01.

Theoretically, the two clades could be classified as two distinct species or two distinct subspecies in D. grignensis. An analysis of the concept and application of the subspecies is beyond the scope of this paper. It is known that subspecies have no place in the schema of the phylogenetic species concept, so adherents to that concept are predisposed to a finding of “no subspecies” because below the generic level a taxon is either a species or it is nothing (Patten 2015). Nevertheless, subspecies are still broadly used, even in papers dealing with cladistic analyses. They also have an importance for conservation issues (Zink et al. 2013), albeit their definition inevitably suffers from a certain arbitrary subjectivity (a charge that could be levied to any rank, of course). Thresholds have been suggested to determine species and subspecies limits with more objectivity (Haig et al. 2008), but in most cases they have been based on a statistical analysis of morphological variation; this has often played a leading role in the definition of a subspecies, indicated as “heritable geographic variation in phenotype.” (Patten 2015). However, morphology does not always seem to be a valid criterion in evaluating taxonomic differences, as demonstrated by the large number of cryptic species that mainly differ in non-morphological traits (Korshunova et al. 2019).

In our opinion, a comprehensive approach to establish the rank to be assigned to allopatric populations should evaluate equally the entire set of morphological, biogeographical and molecular data.

The two clades of the D. grignensis complex differ by several molecular synapomorphies. Any approach to delimit a threshold in interspecific differentiation in genetic divergence is somewhat arbitrary, since it can vary to a great extent among different taxa. A limit of 5% was roughly considered to be a good indicator (Hebert et al. 2003), and a limit of 3% is relatively well working to separate species that are the result of geographic speciation events driven by gradual accumulation of diverging mutations (Sbordoni 2010). This is the case of the two taxa of Dichotrachelus of the D. grignensis complex. Several cases of even lower variation between allopatric, or marginally sympatric yet genetically separated, cryptic species are known, among which the Southern festoon butterflies [Zerynthia polyxena (Denis & Schiffermüller, 1826) and Z. cassandra Geyer, 1828], whose genetic divergence for COI is 1.4% (Zinetti et al. 2013). In previous studies on Dichotrachelus, intraspecific divergence was usually lower than 1%, and reached a maximum of 2.8% for specimens of D. maculosus at the extremes of its range (Meregalli et al. 2018).

Moreover, not all the nucleotide differences between the two groups of populations of the D. grignensis complex are synonymous. The two regions of the COI that were amplified, out of a total of 494 amino acids, differ for 8 between the two taxa (variation of 0.016%). The intraspecific amino acid variation usually peaks near zero (i.e., median divergence is 0.009 for Drosophila flies, Shih et al. 2015), and in the populations at the extremes of the range of D. maculosus Fairmaire, 1869 (Switzerland and Maritime Alps, reciprocally distant several hundred km) only 1 amino acid differs out of 275 in the 825 bp fragment (Meregalli et al. 2018, sequences re-examined). The variation in the COI sequences of the two clades of the D. grignensis complex is relevant, and is confirmed by the variation that is present in the ITS2.

The molecular clock places the disjunction between the demes from the Grigna mountains and those from the Bergamasque Prealps at about 5 Mya. Time for speciation cannot be standardized in any way, since it depends on a number of geographical, ecological and genetic factors (Raboski 2011, and references therein). Several approaches evaluated number of generations and genetic distance required for speciation, discussed by Yamaguchi and Iwasa (2013), even though in most cases speciation was associated with reproductive isolation, which is a simplification of the process (Hausdorf 2011). However, even considering that these species of Dichotrachelus may have a biennial generation time, the more than two million generations in total isolation – a time span determined by the geological events that occurred in the area (see below, chapter 4.2) – would be sufficient for achieving a genetic distance compatible with speciation (Yamaguchi and Iwasa 2013).

We could not evidence any clear-cut morphological trait to separate the two taxa, excluding a minute variation in the shape of the genital sclerite; nevertheless, these minimal differences give an indication that supports the molecular data.

Based on this multi-data approach, we consider the two taxa to be distinct species, even though they are barely distinguishable based on morphology only.

Two main scenarios can be considered. One, the ancestral species was distributed across the entire area from Mount Camino to the Grigna massifs, and then, around 5 Mya, the Grigna population remained isolated from the others. Alternatively, a second scenario supposes that the ancestral species was distributed only in the Bergamasque area, and then, around 5 Mya, it colonized the Grigna massifs, where it remained isolated from the parent population. Some remarks on orogeny and paleoclimatology of the area may help understanding the distribution pattern of the D. grignensis complex. Bioecology of these weevils must be also taken into account. The Dichotrachelus of the D. grignensis group are strictly associated with calcareous rocks; they are, nowadays, apparently exclusive to heights above 1500 m a.s.l., where they are monophagous on Saxifraga, usually S. caesia. Their vagility is very limited and consequently their dispersion capability is scarce, at least in the short term, and it can occur only in suitable habitats with limestone rocks and the presence of Saxifraga caesia. Their present distribution in the Alps was determined by orogeny of the chain as well as paleoclimatic events.

The southern Alps, separated from the Alpine axial zone by the Insubric line, experienced several deformation phases mainly of Late Eocene-Oligocene age (Laubscher 1985; Bersezio et al. 1993). In the Lago di Lecco area, the Grigna Mountains are characterized by thrusted Triassic carbonates (Laubscher 1985). A pre-existing transfer fault zone separates the Grigna block from the Valtorta section, passing through Ballabio and the upper Valsassina (Laubscher 1985).

The southward propagation of thrusting towards the foreland, as recorded in the Po Plain subsurface, continued to occur throughout the Miocene (5.3–7.2 Mya) (Fantoni et al. 2001; Fantoni et al. 2004). Seismic and well data in the Po plain show that shelf margins associated with fluvio-deltaic systems were established at least since the middle Miocene, suggesting an important uplift of the chain and the development of a well-established drainage system linked to the uplifted mountains to the north (Dalla et al. 1992). In addition, there is evidence of Messinian incisions for all the southern Alpine lakes, caused by the relative base-level fall of the Messinian Salinity Crisis (Audra et al. 2007; Rossi et al. 2015). These paleovalleys delivered sediments to fluvio-deltaic systems basinward (Rossi et al. 2015; Rossi et al. 2018), and were later filled following the Pliocene transgression event (Hsü et al. 1977; Rossi et al. 2015). Neoalpine compression in the western Southern Alps likely ended in the Messinian to early Pliocene (Fantoni et al. 2004).

The present Lago di Como topography originated from a deeply incised Messinian canyon during the Messinian Salinity Crisis (Bini et al. 1978), so that uplift and valley deepening preceded the onset of global cooling and Pleistocene glaciations (Audra et al. 2007). In the Como branch of the Lago di Como, the Messinian entrenchment followed a previous paleovalley, whereas the Lecco branch is younger and more significantly imprinted by glaciers (Bini et al. 1978). The Valsassina valley was also incised as a Messinian canyon, which separated the Grigna massifs from the Resegone (Monzini 2013).

The divergence between the Grigna and the Bergamasque Prealps populations of the D. grignensis complex can be roughly dated to around 5 Mya, and there does not seem to be any evidence of a sudden possibility of range expansion westbound towards the Grigna in that epoch. The most likely scenario proposes that the ancestral forms were already present in the entire region during the Miocene, before the Messinian, and became isolated following the Messinian incisions of the Valsassina that separated the Grigna from the Resegone and the other calcareous mountains east of the Introbio valley. Consequently, this primary disjunction between the two lineages of D. grignensis-complex appears to be related to tectonic events.

This pattern of disjunction, with sister species present on the two sides of the Valsassina, is not exclusive to the two species of Dichotrachelus but it occurs also in other scarcely vagile steno-endemisms. Examples are Boldoriella manzoniana Monzini, 1995, from the Grigna massif, with its vicariant B. carminatii (Dodero, 1917) from the Bergamasque Prealps (Monzini 2013), and B. grignensis Monzini, 1987, from the Grigna, with B. focarilei (R. Rossi, 1965) on the east of the Valsassina valley. These, and other, cases probably share the same biogeographical origin.

The populations of D. orobicus present on the limestone outcrops between the Valsassina and the Valcamonica, and probably exclusive to the montane and alpine region, show a limited reciprocal molecular differentiation. Their apparent isolation in the calcareous massifs has undoubtedly a more recent origin, and was likely determined by the Pleistocene climatic oscillations, when glaciers acted as barriers that prevented gene exchange among the populations remained isolated in the nunattaker (Fattorini 2004).