## 2 Materials and Methods
The entire gene was sequenced from 10 samples using a Forward and Reverse strategy for the three overlapping fragments to ensure coverage of all the gene lengths and comparing them with already published sequences of _M. brachycercus_ and _M. savii_ (AY513824.1-AY513828.1; EU158776.1-EU158804.1; KT896517.1-KT896526.1) and _Microtus duodenicmostatus_ (de Selys-Longchamps 1839) (KT896526.1) as outgroups. For the remaining 22 samples, sequencing procedures were limited using a Forward and Reverse strategy on the second and third fragments to allow for comparison with already deposited samples (KT896517.1-KT896526.1). Electropherograms were handled and edited using MEGA 11 [31]. The taxonomic identities of the obtained sequences were evaluated using BLAST+ [32]. Phylogenetic relationships were estimated through Maximum Likelihood (ML) and reconstructions; an ML tree was obtained using MEGA 11 [31], supporting nodes by 1000 bootstrap replicates. Initial trees for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the Jukes-Cantor model and then selecting the topology with a superior log likelihood value. A discrete Gamma distribution was used to model evolutionary rate differences among sites (4 categories (+G, parameter = 0.2461)). The trees were drawn to scale with branch lengths measured in the number of substitutions per site. Phylogenetic analyses were performed on a dataset including the entire set of samples genotyped for the present study (Montefalco specimens) together with the _cytb_ sequences of the _Microtus savii_ group available in GenBank by February 2023 (Supplementary Table S1). Sequences that did not show double signals in the electropherogram were tested through the TCS v.1.21 program [33; 34], setting network the connection limit at 95% and gaps as \"missing\". TCS v.1.21 subdivided the identified haplotypes into haplogroups, describing the relationships between haplotypes/haplogroups found and clarifying the significant number of substitutions connecting haplotypes [35]. The network was visualized and plotted with tcsBU [36]. Statistical selection of best-fit models of nucleotide substitution was performed by means of jModelTest [37; 38]. On the basis of these analysis, the Jukes -Cantor model was used to assess the evolutionary history of the _Microtus_ specimens; Maximum Likelihood and Neighbour Joining methods were inferred estimating standard error with a bootstrap procedure (500 replicates). For the Maximum Likelihood method,a discrete Gamma distribution was used to model evolutionary rate differences among sites (G categories = 4). Spatial or demographic expansion was estimated through [PERSON]'s D neutrality test [39], and rate homogeneity for nucleotide substitutions among lineages was tested through MEGA 11. The risk of analysing nuclear insertions of mtDNA (NUMTs) instead of mitochondrial original sequences was taken into account and prevented using bones as the DNA source [40; 41]. By any chance, prior purification and sequencing, each amplicon was run on 2.5% agarose gel, and the entire sample set clearly showed only a single and discrete band with no evidence of smaller amplicons, excluding the possible presence of NUMTs that are usually shorter than the corresponding mitochondrial gene [42]. Furthermore, NUMTs frequently show several amino acid variations [40] absent in translating the obtained sequences.