We described the genetic variation of cytochrome b gene sequences of blind mole rats in Turkey. We examined 47 individuals belonging to nine cytotypes of three superspecies Nannospalax leucodon, N. xanthodon and N. ehrenbergi in the 402bp gene sequence of cytochrome b. Phylogenetic analyses showed that relationships between cytotypes were well supported, but deeper divergence between species showed insignificant relationships. Cytotypes of N. xanthodon with low diploid number of chromosomes from western Turkey formed a monophyletic group distinct from the populations with higher number of chromosomes (2n = 56-60). The monophyly of N. xanthodon was supported with respect to N. leucodon (2n = 56) in the Bayesian and maximum likelihood phylogenies. The divergence between two analyzed cytotypes of N. ehrenbergi (2n = 52, 2n = 56) was 9.4 %, and the Kilis cytotype (2n = 52) appeared as the basal branch of the whole analysed dataset. N. ehrenbergi cytotypes were paraphyletic and they formed unsupported relationships with previously described N. galili (2n = 52), N. golani (2n = 54), N. carmeli (2n = 58) and N. judaei (2n = 60) from Israel. The results of this study showed that the Nannospalax species complex most likely represents more species than currently recognized, especially in N. xanthodon. We suggest that cytotypes of N. xanthodon and N. ehrenbergi from Turkey should be investigated in detail as possible candidates for being separate species.
The mole rats are adapted for subterranean life. They are distributed in the Palaearctic region, throughout eastern and southeastern Europe, Anatolia, the Caucasus, and the Middle East up to northeastern Africa (Topachevskii 1969, Wilson & Reeder 1993). Their evolutionary history and taxonomic status are difficult to ascertain and molecular, karyological and morphological studies are needed to determine the phylogenetic relationships of mole rats in Turkey (Nevo et al. 1995, Suziki et al. 1996, Sözen et al. 2000, Kankılıç et al. 2005, Ivanitskaya et al. 2008). For several decades, scientists agree that taxonomy of mole rats (Spalacinae, Rodentia) needs a modern revision based on chromosome and molecular genetic data coupled with morphology, physiology and behavior (Savić & Nevo 1990, Kryštufek & Vohralík 2009). The genus Spalax differs from representatives of the genus Nannospalax in having a less variable diploid number (2n = 60, 62 in Spalax versus 2n = 36-60 in Nannospalax) and a higher number of subtelocentric chromosomes (NF = 116-124 in Spalax versus NF = 72-98 in Nannospalax cytotypes) (Lyapunova et al. 1971, Savić & Nevo 1990, Németh et al. 2009), and shows slow chromosomal evolution rate (Kryštufek & Vohralík 2009). Kryštufek & Vohralík (2009) and Németh et al. (2009) ranked Nannospalax species as superspecies. Here we use superspecies order for Turkish mole rats instead of species.
According to Wilson & Reeder (1993) and Yiğit et al. (2006), three species, Nannospalax leucodon, N. xanthodon (senior synonym of N. nehringi; Kryštufek & Vohralík 2009) and N. ehrenbergi, occur in Turkey. N. leucodon is found in the Turkish Thrace, and N. ehrenbergi in southeastern Turkey. N. xanthodon has a wide distribution area extending over Anatolia.
Karyological studies revealed 17 cytotypes in Turkey: one cytotype (2n = 56) in N. leucodon, 11 (2n = 36, 38, 40, 48, 50, 52, 54, 56, 58, 60 and 60R) in N. xanthodon and five (2n = 48, 52, 54, 56, 58) in N. ehrenbergi (see reviews in Sözen et al. 1999, 2006a, 2011, Coşkun et al. 2006, Kankılıç et al. 2007, 2010). Additionally, Nevo et al. (1994, 1995) recorded the 2n = 62 form, but Ivanitskaya et al. (2008) reported that the 2n = 62 forms should be eliminated from the list of Turkish mole rats. Later Arslan et al. (2011) studied C- and AgNOR banding patterns of three cytotypes (2n = 40, 58 and 60) from southern Anatolia. These karyological results show that one of the most complex chromosomal diversity within the distribution range of the genus Nannospalax is found in Turkey. These results complicate the taxonomical status of the cytotypes in Nannospalax in Turkey.
While taxonomic status of cytotypes of mole rats in Turkey has been under discussion, Nevo et al. (2001) raised four cytoypes (2n = 52, 54, 58 and 60) in Israel to species level; Nannospalax galili (2n = 52), N. golani (2n = 54), N. carmeli (2n = 58) and N. judaei (2n = 60). Arslan et al. (2010) studied mitochondrial DNA (mtDNA) variation between three cytotypes (2n = 40, 58 and 60) in southern Anatolia, and they showed well-supported lineages in the phylogenetic tree. But there is no detailed study based on mtDNA sequences on the genetic structure of cytotypes of mole rats in Western Turkey. This study focused on phylogenetic relationships among cytotypes of N. xanthodon in Western Turkey, and N. leucodon in the Turkish Thrace and N. ehrenbergi in southeast of Turkey inferred from mtDNA cytochrome b gene sequences. We aimed to highlight the genetic relationships between and among the cytotypes of the species in Turkey, and also the recently described species of Nannospalax from Israel.
Material and Methods
Sampling for molecular studies
A total of 47 mole rat individuals were used in the molecular studies. These samples came from three previously recognized species, namely N. leucodon (N = 3), N. xanthodon (N = 40; cytotypes 2n = 36, 38, 40, 50, 56, 60) and N. ehrenbergi (N = 4; cytotypes 2n = 52 and 56) (Fig. 1). Karyotypes were prepared according to Ford & Hamerton (1956). Skins and skulls were deposited to the Department of Biology of Zonguldak Karaelmas and Ankara University, Turkey.
DNA isolation, amplification, and sequencing
We sequenced all 47 samples for a part of the cytochrome b locus. Total genomic DNA was isolated from muscle tissues following the techniques reported by Doyle & Doyle (1987) known as CTAB method. Fragment of the cytochrome b gene was amplified using polymerase chain reaction in 25 µl reaction volume and each reaction included 1 µl of each primer (20 pmol) (L14724 and H15154; Irwin et al. 1991, Smith & Patton 1993), 4 µl of dNTPs, 2.5 µl of 10× Taq Buffer with (NH4)2SO4, 1.5 µl of MgCl2 (25 mM) and 0.25 µl of Taq DNA Polymerase (5 units/µl, Fermentas, Ontario, Canada). The PCR cycling conditions first included four cycles of 1 min of denaturation at 95 °C, 1 min of annealing at 40 °C, and 1 min extension at 72 °C, followed by 33 cycles using annealing temperature at 50 °C. Before the sequencing reaction the PCR fragments were cleaned with Nucleospin extract kit (Macherey-Nagel, Düren, Germany) and later sequencing reactions of both DNA strands were commercially performed using Big Dye Terminator v. 3.1 sequencing chemistry on an ABI 3100 Genetic Analyzer (Applied Biosystems, California, USA). All specimens were deposited to the GenBank under the accession numbers FJ656259-FJ656305 (Table 1).
Samples used in the present study.
Previously published cytochrome b sequences of Nannospalax were obtained from GenBank (AF140523-AF140534; Nevo et al. 1999), and two other sequences for Mus musculus (NC010339) and Rattus norvegicus (EU273707) were added as an outgroup. The dataset was aligned in Clustal X 2.0.9 (Larkin et al. 2007) and had the total length 402 basepairs.
Genetic divergence was estimated in MEGA 4 (Tamura et al. 2007). Neighbour-joining (NJ) and maximum parsimony (MP) phylogenetic analyses were executed in PAUP 4.0b10 (Swofford 1999), using HKY + G substitution model for the NJ analysis selected by Bayesian Information Criterion in Modeltest 3.7 (Posada & Crandall 1998), which was the simplest suitable model, and 10000 bootstrap replicates. Maximum likelihood (ML) tree was obtained from RAxML 7.4 (Stamatakis 2006) with GTR + G model and 10000 bootstrap replicates. Bayesian Inference (BI) analysis was run in MrBayes 3.1.2 (Ronquist & Huelsenbeck 2003) with 12 Markov chains Monte Carlo for two million generations in two independent runs. The chain swapping was successful in 68-72 % of attempts indicating effective chain mixing. We used default prior settings with the exception of the branch length. That was limited to brlenspr = unconstrained: exp(20) to ensure that the posterior 95 % confidence interval of the tree length included the tree length obtained from the maximum likelihood analysis (cf. Marshall 2010). The prior setting did not influence tree topology compared to the default setting (data not shown). A median-joining network (Bandelt et al. 1999) was also reconstructed for the sequenced mtDNA region. Finally, tree topologies were compared by using the Shimodaira-Hasegawa test in PAUP 4.0b10 (Swofford 1999).
The 47 sequences represented 42 unique haplotypes, which were used in the phylogenetic analyses. The 402 bp long dataset contained 120 parsimony informative sites and 31 singleton mutations. Support for monophyly of the studied taxa based on partial cytochrome b sequences was different for specific analyses (Fig. 2). Total tree length in the MP analysis was 364 steps and the groupings reflected the genetic distances between cytotypes and species.
Nucleotide diversity within cytotypes ranged from 0 to 0.03 substitutions per site (Table 2). Within N. xanthodon, the greatest diversity was found in the 2n = 38 cytotype that was also represented by the highest number of individuals in our dataset.
The genetic distances were calculated between all taxa used in the present study using uncorrected p-distance (Table 3). The distance between four Israeli Nannospalax species ranged between 0.019 and 0.058 substitutions per site. The cytotypes attributed to N. xanthodon showed deeper divergence ranging between 0.024 and 0.105. The greatest genetic distance between studied cytotypes was between N. carmeli and N. xanthodon 2n = 40. Smaller genetic distances tended to be found between the chromosomal cytotypes with the highest diploid chromosome numbers. The genetic distance between the two cytotypes of N. ehrenbergi in Turkey was 0.087, which is greater than the genetic distance between the two most distantly related Israeli species (N. golani and N. carmeli) (Table 3).
Within group nucleotide diversity for partial sequences of the mitochondrial cytochrome b gene (402 bp) for identified Nannospalax cytotypes.
Pairwise uncorrected p-distance for partial sequences of the mitochondrial cytochrome b gene (402 bp) for identified Nannospalax cytotypes (below diagonal) and their standard errors (above diagonal).
We conducted four separate sets of analyses. The data were analyzed using Bayesian inference, neighbourjoining, maximum likelihood and maximum parsimony analyses. The trees showed similar topologies (Fig. 2). The result of the Shimodaira-Hasegawa test showed that only NJ tree topology was significantly different from the other trees (Table 4). This result was caused by the placement of the N. ehrenbergi clade. In all trees, cytotypes with low diploid chromosome numbers (2n ≤ 50) grouped as one supported clade.
Results of Shimodaira-Hasegawa (SH) tests for all trees. The NJ tree topology was found different from ML (* p < 0.05).
The other N. xanthodon cytotypes (2n = 56 and 60) occurred as a separate clade, but monophyly of N. xanthodon was supported only in BI and ML analyses (Fig. 2). N. leucodon formed a monophyletic group with N. xanthodon. The other taxa showed a basal polytomy in most of our analyses. Both cytotypes of N. ehrenbergi were paraphyletic with respect to the Israeli taxa. N. carmeli and N. judaei formed a monophyletic group, but the taxa were not reciprocally monophyletic within this group.
We found similar topology in the network analysis (Fig. 3). The N. xanthodon cytotypes formed two putative groups with low (2n ≤ 50) and high (2n ≥ 56) chromosome numbers. The former group was more closely related to N. leucodon as well as to the group including N. ehrenbergi and the Israeli taxa (Fig. 3).
Our phylogenetic analysis showed monophyly of N. xanthodon and N. leucodon, but paraphyletic relationships of the two N. ehrenbergi sequences with Israeli taxa. N. xanthodon sequences formed two reciprocally monophyletic groups that included cytotypes with low (2n ≤ 50) and high (2n ≥ 56) chromosome numbers. Their sister species was N. leucodon. These results are similar to recent findings of Arslan et al. (2010) on a smaller dataset including three cytotypes and Kryštufek et al. (2012) that were based on a longer alignment. The taxa N. ehrenbergi, N. galili, N. golani, N. carmeli and N. judaei diverged in a basal polytomy wherein only a single group including N. carmeli and N. judaei was consistently retrieved. While Reyes et al. (2003) also showed an undifferentiated clade of N. carmeli and N. judaei samples, our analyses did not support their clear relationship of N. galili and N. golani. This was distorted by position of N. ehrenbergi 2n = 56 as a sister taxon with unresolved relationship to the Israeli species.
Cytotypes of N. xanthodon with low diploid number (2n = 36, 38, 40 and 50) formed a monophyletic group. This may be a result of specific evolutionary pathways. For example, Matur & Sözen (2005) stated that the River Sakarya separated 2n = 52 and 60 cytotypes of N. xanthodon in Bilecik province, and the river might act as a barrier between these cytotypes. But in our analyses, 2n = 52 N. xanthodon was not analysed and we cannot confirm this scenario. According to our results, 2n = 38 from Gökçeada (Aegean island) differentiated from the 2n = 38 populations in the mainland indicating mole rat diversification in island isolation.
Monophyly of specific cytotypes was present in all cytotypes belonging to the N. xanthodon 2n ≤ 50 group, but the cytotypes were not likewise differentiated in the 2n ≥ 56 group. Rather, the northern populations of N. xanthodon, Karabük (2n = 56N) and Kastamonu (2n = 60R), were differentiated from the other samples of 2n = 56 and 60. This was suggested by Sözen (2004) and Sözen et al. (2006b), who claimed that the cytotypes from Northern Anatolia (2n = 54N, 56N, 58N and 60R) differentiated from other cytotypes. Ivanitskaya et al. (2008) showed that the 2n = 60 population in Kastamonu is different from the 2n = 60 population in central Anatolia using classical and molecular cytogenetic techniques. Matur et al. (2010) studied chromosomal evolution of Turkish mole rats inferred from G- and C- banding of 11 cytotypes. They found that cytotypes from Kastamonu (2n = 60R) and Karabük (2n = 56R) formed a clade together with 2n = 54 cytotype not included in our study. This northern clade has independent evolutionary pathway, and 2n = 60R might be considered an ancestral karyotype (Matur et al. 2010).
There is a conflict in the ancestral population of Thracian mole rats (N. leucodon 2n = 56). They grouped with N. xanthodon in our trees. Recently, Matur et al. (2011) proposed the 2n = 60 cytotype as an ancestral for all Anatolian cytotypes. Nevertheless, this process needs more detailed studies in order to verify the ancestral karyotype of the N. xanthodon/leucodon group.
In our phylogenetic analyses, N. carmeli (2n = 58) and N. judaei (2n = 60) did not appear in separate clusters. Reyes et al. (2003) found similar results from the sequence studies of the mitochondrial control region. In their study, N. golani and N. galili were well separated, whereas N. judaei (2n = 60) and N. carmeli (2n = 58) formed a single cluster. Nevo et al. (1999) argued that this is due to the fact that they represent young species. The paraphyly of the two cytotypes of N. ehrenbergi (2n = 52 and 2n = 56) from Kilis and Osmaniye further complicate this in our study. Both cytotypes formed unsupported relationships and polychotomies. Tentatively, N. ehrenbergi from Osmaniye grouped with N. galili and N. golani, and N. ehrenbergi from Kilis was a basal taxon in the tree (Fig. 2).
Further comprehensive and detailed multidisciplinary research combining morphology, karyology, physiology, behaviour, mtDNA and nuclear DNA phylogeny should be applied to clarify the taxonomic status and phylogenetic relationships of cytotypes of N. leucodon, N. xanthodon, and N. ehrenbergi in Turkey. In addition, research on intra-population variation should also be considered for better understanding of mole rats in Turkey.
This study was supported by the Scientific and Technical Research Council of Turkey (TUBITAK 101T084 and 106T225), Zonguldak Karaelmas University (2004-13-06-08, 2008-13-06-01) and Academy of Sciences of the Czech Republic (AV0Z60930519). The authors thank Ahmet Turkyılmaz from REFGEN Biotechnology Inc. for performing the sequencing. Authors thank to Zeynep Yuksel Sezen and Joan E. Johnston for careful editing the manuscript and two anonymous reviewers and the editor for helpful and detailed comments.