Samples

This study consists of seven snow voles collected in three regions of Iran (Table 1 and Fig. 1), two of which were topotypical with C. layi. All specimens are preserved as standard museum skins and skulls and deposited at the Museum of Vertebrate Zoology, Berkeley, U.S.A. (MVZ) and the Zoological Museum of Ferdowsi University of Mashhad, Iran (ZMFUM, Table 1). Taxonomic identification was facilitated by examination of traditional taxonomic characteristics (colour, external and cranial measurements, and enamel molar pattern) provided by Kryštufek & Vohralík (2005) and Zykov (2004). Cranial measurements were scored by Vernier calliper to the nearest 0.1 mm and the external measurements were deduced from specimen tags.

Fig. 1.

Records of snow voles Chionomys in Iran: 1 — Central Alborz Mts., 2 — Binaloud Mts., 3 — Zagros Mts. See Table 1 for further locality information. Range (shaded grey) is modified from Kryštufek & Amori (2008).

Cytb sequences of our voles were combined with 94 sequences downloaded from GenBank (supplementary information, Table S1) and belonging to all three currently recognized species of snow voles (Jaarola et al. 2004, Pfunder et al. 2004, Galewski et al. 2006, Bužan & Kryštufek 2008, Castiglia et al. 2009, Bannikova et al. 2010, 2013, Fink et al. 2010, Yannic et al. 2012, Barbosa et al. 2013, Arslan et al. 2016). Phylogenetic trees were rooted using three arvicolines which are closely related to Chionomys (Zou et al. 2008, Bannikova et al. 2009): Blanfordimys afghanus (GenBank Accession Numbers EF599108-9), B. juldaschi (EF599112-3), and Microtus fortis (EU126809). To exclude possibility of introgression of alien mtDNA via hybridisation (Harrison & Larson 2014), we sequenced the nuclear gene for growth hormone receptor (ghr) in two topotypes of C. layi and two individuals from Binaloud Mts. All available ghr sequences from GenBank were used as supplementary data (C. nivalis: AM392378, JN244718-9, JX440339, JX440340, KP057338_ C. gud: JN244708-14_ C. roberti: JN244715-17; Galewski et al. 2006, Bannikova et al. 2013, Petrova et al. 2016). Arvicola amphibius (ZMFUM_4397), Blanfordimys afghanus (ZMFUM_1426), Microtus pennsylavnicus (AF540633, Adkins et al. 2003) and Ellobius talpinus (PMS_223) were used as putative close relations to rooting phylogenetic trees.

Table 1.

Details of sample localities (mapped in Fig. 1) for snow voles from Iran. Geographic origin of the samples and GenBank accession numbers are provided. Museum vouchers are deposited at the Museum of Vertebrate Zoology, Berkeley, U.S.A. (indicated by asterisk *) and the Zoological Museum of Ferdowsi University of Mashhad, Iran (the rest).

Molecular analyses

Nucleotide sequences were edited manually in CodonCode Aligner software (CodonCode Corp.), and aligned with the Clustal W (Thompson et al. 1997) algorithm using BioEdit 7.0.5 (Hall 1999). The best fit model of DNA substitution was chosen using jModeltest 0.1.1 (Posada 2008), based on the Akaike Information Criterion (AIC). The General Time-Reversible (GTR) model with proportion of invariant sites (I = 0.5502) and gamma distribution of rates across sites (G = 1.1005) best fitted cytb dataset. The best fit model for ghr dataset was HKY + G (G = 0.2193). Pairwise genetic distances, both between species and within them, were estimated with Kimura-2 parameter model (K2P) using Mega v6 (Tamura et al. 2013).

Fig. 2.

(A) Bayesian inference tree reconstructed from cytochrome b sequences of the genus Chionomys and rooted on Microtus fortis, Blanfordimys afghanus and B. juldaschi. (B) BI tree generated from the ghr sequences of the genus Chionomys. The numbers on the branches correspond to posterior probability values (BPP > 0.90) and bootstrap supports (BP > 70 %). The triangles represent species and lineages of snow voles which are based entirely on published sequences. For haplotype designations see Table 1 and supplementary information (Table S1).

Table 2.

K2P genetic distances between the geographical lineages of the European snow vole Chionomys nivalis. Within lineage distances are shown in diagonal.

Bayesian inference (BI) and Maximum Likelihood (ML) phylogenetic trees were constructed using selected substitution model. BI analysis was conducted in MrBayes 3.1.2 (Huelsenbeck & Ronquist 2001, Ronquist & Huelsenbeck 2003), using two simultaneous analyses with four Markov Chain Monte Carlo (MCMC) models which started from random trees and were run for 4 million iterations. The trees were sampled every 1000th generation after removing the first 10 % of the trees as the burn-in stage. ML tree estimation was carried out using PAUP 4.0b10 (Swofford 2003) under 100 bootstrap pseudo-replicates with a ML heuristic tree search using 10 random additional sequence replicates. Nodal robustness for BI and ML analyses were assessed using Bayesian Posterior Probability (BPP) and Bootstrap values (BP) for BI and ML analyses, respectively. The values BPP > 0.95 and BP > 80 were assumed as good support for BI and ML analyses, respectively.

Taxonomic identity

Comparison with published character states classified our specimens as snow voles (Chionomys). Any of them showed similarities with the Caucasian endemics C. gud and C. roberti but matched C. nivalis-layi. Zykov (2004) distinguished C. layi from C. nivalis by longer ears (20 mm), longer diastema (above 9.5 mm) and wider bullae (breadth over 7 mm). Morphological comparison with the data provided by Zykov (2004) showed broad similarity between the topotypes and the type series of C. layi. Corresponding measurements in our topotypes of layi were (parenthesized is range for C. nivalis in Zykov 2004): length of ear 19 and 20 mm (12-18 mm), length of diastema 9.0 and 9.7 mm (7.5-9.3 mm), and breadth of bulla 7.0 and 7.4 mm (6.0-7.0 mm).

Molecular analyses

Altogether, five new cytb haplotypes were found in our material generating 74 different snow vole cytb haplotypes. Within the 1140 bp long sequences considered here, 193 polymorphic sites were found with a total of 267 mutations, 170 of which were parsimony informative. Our analysis of ghr retrieved three new sequences (one from Zard Kuh and two from Binaloud Mts.) generating 20 different ghr sequences. Within 870 bp long sequences, 25 and 20 sites were found to be polymorphic and parsimony informative, respectively. BI and ML trees of mitochondrial data showed similar topologies, therefore only the BI tree is shown (Fig. 2A). Branching pattern in both trees was congruent with published results (Bannikova et al. 2013, Arslan et al. 2016) therefore three species of snow voles emerged also in our reconstruction: C. nivalis, C. gud, and C. roberti. Sister position of C. nivalis against the clade of C. gud + C. roberti benefited strong BP support (100 %). Within C. nivalis the enigmatic sequence from Aladağ (JN244707, Bannikova et al. 2013) hold a supported (BP 95 %) basal position. The remaining haplotypes clustered into two strongly supported (BPP 1.00, BP 99 %) major lineages which contained also our new sequences from Iran. Four samples from Binaloud Mts. hold a supported (BPP 1.00, BP 99 %) sister position against the reference samples of C. n. dementievi from adjacent Kopet Dagh in Turkmenistan (Turkmenian lineage). The remaining haplotypes from Europe and Asia were in a single cluster (Euro-Asiatic lineage). Firstly, two strongly supported (BPP 1.00, BP 89 %) sublineages emerged in this lineage, the European and the Asiatic. Further substructuring within the latter retrieved two supported (BPP 0.98) clusters, the Caucasian and the Anatolian. The latter contained published samples from Turkey and the Levant coast in Syria and Israel. Our samples from Alborz (= Elburz) and Zagros clustered with the Caucasian haplotypes. This cluster however further substructured (BPP 0.98) and the two sequences from the Zagros were in a sister position against samples from Alborz + Caucasus. Phylogenetic analysis of ghr gene showed identical tree topology (Fig. 2B).

Putting aside the enigmatic haplotype from Aladağ, Turkey, K2P distances between European snow vole lineages varied between 1.8 ± 0.4 (Syria + Israel vs. Anatolia) and 5.7 ± 0.7-6.1 ± 0.9 (Turkmenian vs. the remaining lineages); median K2P was 3.8 %. Within-lineage distances were mainly below 1 %; higher values (1.3 ± 0.3-1.7 ± 0.3) evidently point on further substructuring, particularly so in the European lineages, and within the Caucasian lineage. Samples from Zard Kuh were separated from a cluster Alborz + Caucasus by K2P distance of 3.2 ± 0.5 (Table 2).