Phylogenetic Status of the Turkish Red Fox ( Vulpes vulpes ) , based on Partial Sequences of the Mitochondrial Cytochrome b Gene

Genetic diversity and multiple mitochondrial phylogroups of the red fox have been revealed from scattered locations in previous studies. There is a still lack of information about the genetic diversity and phylogeographic structure of the red fox in Asia Minor. We investigated the genetic diversity in the Turkish red fox using a part of the cytochrome b mitochondrial gene (375 bp), and attempted to evaluate the phylogeographic structure in various geographic ranges of the species with the use of sequences available from the GenBank from various geographic origins and our data. Bayesian and Network analyses of the cytochrome b sequences from Turkey and GenBank suggested that the red fox is divided into four main phylogroups. They are grouped accordingly: Group 1 (SW Anatolia, Turkey and Hokkaido, Japan), Group 2 (Eurasia and North America), Group 3 (only North America), and Group 4 (Vietnam). The majority of Turkish haplotypes grouped with those of Eurasia. Despite the great distance between the localities, two haplotypes from SW Anatolia, Turkey grouped with previously reported haplotypes from Hokkaido, Japan. The present study shows that the Turkish red fox is nested within two main phylogroups and exhibits a high genetic diversity.

Multiple studies based on mitochondrial DNA of red fox emphasize that there has been a considerable importance to have a comprehensive picture of genetic variation in the entire range of the species, when new insights are added into phylogeographic and colonization patterns of red fox with a wider biogeographic context (Frati et al. 1998, Inoue et al. 2007, Perrine et al. 2007, Aubry et al. 2009, Teacher et al. 2011, Edwards et al. 2012, Yu et al. 2012a, Kutschera et al. 2013).Some of the studies resulted in the formation of mitochondrial phylogroups, which were classified as Holoarctic and/or Eurasian, Nearctic, and Hokkaido (Inoue et al. 2007, Aubry et al. 2009, Yu et al. 2012a, Kutschera et al. 2013).Although there is a substantial inventory of mitochondrial data from Eurasia and North America (Frati et al. 1998, Inoue et al. 2007, Perrine et al. 2007, Aubry et al. 2009, Teacher et al. 2011, Edwards et al. 2012, Yu et al. 2012, Kutschera et al. 2013), there is a still lack of information about the genetic diversity and phylogeographic structure of the red fox from Turkey and the Middle East, except Israel.In this sense, more genetic and phylogeographic studies are needed to place these new mitochondrial data in a wider context and to elucidate the current pattern of phylogeographic structure throughout the red fox' range.Therefore, it would not only be beneficial to determine the genetic structure of different red fox populations but also to focus on regional details, especially those in previously unsampled areas.Turkey is one of such areas.Because of its continental position, Turkey constitutes a mysterious and interesting area for mammalian biogeography due to its climate, vegetation, heterogeneous topography and geological history; together these unique factors produce a biodiversity hotspot (médail & Quézel 1999, Bilgin 2011).Currently, there are at least 150 or more mammalian species from Turkey (Kryštufek & Vohralik 2001, 2009, Wilson & Reeder 2005).Although there are a few studies based on its morphology and distribution (Özkurt et al. 1998(Özkurt et al. , Temizer 2001)), the phylogenetic status and genetic diversity of Turkish red fox have not been previously investigated.Therefore, there is a lack of information in this respect.In reference to the red fox dentition study, Szuma (2008) suggested that the modern red fox evolved either in Asia Minor (Ana tolia) or North Africa during the early Pleistocene.Con sequently, to confirm this, the phylogenetic status of Tur kish red fox should be revealed using different markers.
In this study, our primary purpose was to determine the magnitude of genetic diversity within the Turkish red fox by sequencing part of the Cyt b gene.With emphasis on this information we provide new insights into phylogeographic processes which may be essential for a better understanding of evolution.

Localization, tissue sampling and DNA extraction
Various tissues (ear, tail, skin and skeletal muscle) of 51 red foxes, which were road kills, were collected from 51 different localities in Turkey (Table 1, Fig. 1) during 2002-2010.Tissue samples had been stored in the corresponding author's private collection at the Department of Biology, Faculty of Sciences, Erciyes University in Kayseri, Turkey.
Total genomic DNA was extracted by using the DNeasy Blood and Tissue Kit (QIAGEN), following the manufacturer's instructions.A segment of the mi tochondrial cytochrome (Cyt) b (375 base pairs) from each Turkish sample was amplified by polymerase chain reaction (PCR) using the primer pair L14724 (5' GATATGAAAAACCATCGTTG 3') and H15149 (5' CAGAATGATATTTGTCCTCA 3')  (see Kocher et al. 1989, Irwin et al. 1991, Frati et al. 1998).PCR amplifications were performed in 50 µl reaction mixture (1 × Taq buffer with (NH 4 ) 2 SO 4 , 200 µM dNTP mix, 2 u Taq DNA polymerase (Fermentas), 2 mM MgCl 2 , 1 µM of each primer, 1 µl DNA extract).The PCR condition was 3 min at 94°C, 30 cycles of 1 min at 94°C, 1 min at 50°C, and 1.5 min at 72°C, and 7 min at 72°C.The quality of the total DNA and the PCR products were verified by electrophoresis in 0.8% and 2% agarose gels, respectively, and stained with ethidium bromide.PCR products were purified using High Pure PCR Product Purification Kit (Roche).Sequencing was performed in both directions (Forward and Reverse) using the PCR primers, resulting in sequences readable for 375 bp.The purified PCR products were sequenced with an ABI 3100 Genetic Analyzer (RefGen, METU, Technopark-Ankara, Turkey).

Phylogeographic analyses of mtDNA haplotypes of red fox
To place these new mitochondrial DNA data in a wider context, we included all previously published red fox Cyt b gene sequences in the phylogenetic analyses.To the best of our knowledge, this represented all matrilineal phylogeographic studies from various parts of red fox range in the North Hemisphere (Ledje & Arnason 1996 Teacher et al. (2011) and Statham et al. (2012) were slightly shorter (354 bp and 250 bp).Therefore, the sequences were truncated to align to the shorter published sequences (354 bp and 250 bp).
For all datasets, phylogenetic relationships between haplotypes of the Cyt b gene were revealed by using median-joining (MJ) network within the software Network 4.6.1.1 (Bandelt et al. 1999, http://www.fluxus-engineering.com) to represent the intra-specific genealogy of the haplotype datasets.To solidify statistical analysis for the mitochondrial Cyt b phylogroups, we also constructed phylogenetic trees with Bayesian phylogenetic analysis (BI: Bayesian inference of phylogeny) as implemented in the program Mr.Bayes v.3.2.1 (RonQuist & Huelsenbeck 2003).For all datasets, the HKY (Hasegawa-Kishino-Yano) + G (gamma) model of nucleotide substitution was determined as the most appropriate model according to Bayesian Information Criterion (BIC) using the program jModeltest v.0.1.1 (Posada 2008).The four Monte Carlo Markov chains were used to calculate the Bayesian posterior probabilities for 2 million (for 250 bp), 4 million (for 354 bp) and 2 million (for 375 bp) generations with trees sampled every 100 generations; the first 25 % of samples were discarded as burn-in (average standard deviation of split frequencies <0,01).After discarding burn-in trees and evaluating convergence, the remaining samples were retained for generating consensus trees (50 % majority rule), calculating 95 % Bayesian credibility intervals and posterior probabilities.Bayesian tree diagrams were obtained with tree figure drawing tool, FigTree v1.3.1 (Rambaut 2009).Additionally, we obtained genetic divergence among the phylogroups by using K2P (Kimura's 2-parameter distance) estimator (Kimura 1980) implemented in MEGA5 (Tamura et al. 2011).An approximate estimation of divergence time among the phylogroups was calculated on the basis of genetic divergence obtained with K2P, using Avise (2000)'s formula: where P net is the corrected distance between the A and B phylogroups, P AB is the mean genetic divergence in pairwise comparisons of haplotypes A versus B, and P A and P B are the mean genetic divergence among haplotypes within these phylogroups.We also used the Cyt b divergence rate of 2 % per million years as large mammals (Brown et al. 1979, Avise et al. 1998).The arctic fox, Alopex lagopus, [Accession number: AY598511 (Delisle & Strobeck 2005)] was used as out-group following Teacher et al. (2011).

Phylogenetic relationships among the red fox haplotypes
We identified a total of 54, 79 and 45 haplotypes in the three data sets: 250 bp, 354 bp and 375 bp the Cyt b sequences, respectively, including sequences from Gen-Bank and Turkey (51 samples).The three data sets were analyzed by using the Bayesian phylogenetic analysis as the probabilistic method and the MJ network analysis.The two datasets, 250 bp and 354 bp, displayed similar topology in trees and networks and, therefore, here we presented only the Bayesian tree and MJ network for 250 bp sequences.
Based on the both 250 bp and 354 bp Cyt b sequences, the MJ networks formed four distinct phylogroups (Fig. 2): Group 1 (Southwest Anatolia, Turkey and Hokkaido, Japan), Group 2 (Eurasia and North America), Group 3 (only North America) and Group 4 (Vietnam).Based on the 375 bp Cyt b sequences from Eurasia and Turkey, the MJ network formed in two main distinct phylogroups (Fig. 3): Group 1 (The SW Anatolia, Turkey and Hokkaido, Japan) and Group 2 (Eurasia).
In addition to the MJ network analysis, we also used Bayesian phylogenetic analyses for the three data sets (250 bp, 354 bp and 375 bp).Bayesian analyses revealed that the general picture of red fox phylogeography was in agreement with those of previous studies.
In the two data sets (250 bp and 354 bp) that were analyzed, the majority of red fox haplotypes were divided into three major phylogroups as reported by Inoue et al. (2007), Aubry et al. (2009), Yu et al. (2012a) and Kutschera et al. (2013) (Fig. 4).Of the 10 Turkish haplotypes, TR7 and TR10 clustered within the SW Anatolia -Hokkaido phylogroup (Group 1), and the remaining Turkish haplotypes grouped into Group 2 (Eurasia and North America).Furthermore, the Bayesian analysis based on the 375 bp Cyt b haplotypes formed two distinct groups; Group 1 (SW Anatolia, Turkey and Hokkaido, Japan) and Group 2 (Eurasia) (Fig. 5).In all Bayesian trees (Figs. 4, 5), the phylogeographic patterning of haplotypes in the three data sets exhibited that   , 5), the distribution pattern of red fox haplotypes did not exhibit a strong geographical dispersion at continental level.By using the mean genetic distances calculated by K2P estimator (Kimura 1980) and by taking the suggested Cyt b divergence rate of 2% per million years as large mammals (Brown et al. 1979, Avise et al. 1998), we estimated that the divergence time on the basis of the 250 bp Cyt b sequences was around 1.68 million years ago (mya) between Group 1 and 2, and 0.99 mya between Group 1 and 3 (Figs.2, 4), on the basis of the 354 bp Cyt b sequences was around 1.27 mya between Group 1 and 2, and 0.71 mya between Group 1 and 3 (not shown Figure ), and on the basis of the 375 bp cyt b sequences was around 1.22 mya between Group 1 and 2 (Figs. 3, 5).

Discussion
By using a molecular marker, mitochondrial DNA, and a large number of samples from Turkey, we were able to present the first considerable work to assess the ge-netic diversity and phylogenetic status of the Turkish red fox.Our results (Figs. 2 -5) verified all phylogroups (clades) of red fox described by previous studies (Inoue et al. 2007, Aubry et al. 2009;Yu et al. 2012a): (i) Group 1 was spread within Hokkaido, Japan and SW Anatolia, Turkey; (ii) Group 2 in Eurasia and North America; (iii) Group 3 in North America, and (iv) Group 4 consisted of one haplotype from Vietnam.Our study showed that the Turkish red fox was nested within two main phylogroups (Figs. 2 -5).The divergence time between phylogroups of red fox estimated on the basis of the mitochondrial Cyt b variations were found to be around 0.71 -1.68 mya, and had occurred during different periods of the Pleistocene.Except for one data set (250 bp), compared to previous studies (Inoue et al. 2007, Yu et al. 2012a, Kutschera et al. 2013), our estimations were compatible with those of Inoue et al. (2007) using Cyt b and control region, and Yu et al. (2012a) using Cyt b, but in fact quite different from those of Kutschera et al. (2013), who used the mitochondrial control region.A similar depiction of divergence times was drawn for some large animals e.g., brown bear and red deer as in the studies carried out by Matsuhashi et al. (1999Matsuhashi et al. ( , 2001)), Leonard et al. (2000), Mahmut et al. (2002) and Skog et al. (2009).The above mentioned mammals (Ursus arctos and Cervus elaphus) and red fox (Vulpes vulpes) hereby are assumed to share a similar evolutionary history during the Pleistocene (Yu et al. 2012a).
By comparing the haplotype and nucleotide diversity in Turkish samples from the studied region with data from other parts of the species range, we found ten distinct haplotypes (375 bp) in the 51 Turkish samples on the basis of partial sequences of the mitochondrial Cyt b, (with the average one different haplotype over 51/10=5.1 samples) and haplotype diversity=0.7380.Frati et al. (1998) found 18 distinct haplotypes (375 bp) in 41 samples from the Mediterranean Basin (with the average one different haplotype over 2.28 samples).Although Inoue et al. (2007) found 14 different haplotypes (375 bp) in 88 samples from Northern Japan (with the average one haplotype over 6.28 samples), Yu et al. (2012a) determined nine haplotypes (338 bp) in 22 samples from East Asia (with the average one different haplotype over 2.44 samples).By analyzing the North American red foxes, Perrine et al. (2007) and Aubry et al. (2009) found 14 haplotypes (354 bp) in 73 samples (with the average one different haplotype over 5.21 samples), and 29 haplotypes (354 bp) in 220 samples (with the average one different haplotype over 7.58 samples), respectively.In conclusion, the Turkish red fox exhibits a high genetic diversity.If study areas are compared, the scale of our study area was geographically much smaller than those of Frati et al. (1998), Aubry et al. (2009) and Yu et al. (2012a), but much larger than those of Perrine et al. (2007) andinoue et al. (2007).Of the 10 Turkish haplotypes from the 51 samples using a part of the Cyt b gene (375 bp), TR1, TR2, TR4, TR5, TR6, TR8, TR9 and TR10 were new haplotypes for red fox, and may be local haplotypes found in Turkey.TR1 was a common haplotype in Turkey (Table 1; Fig. 1), with an ecologically monotonous (step area and plateau) region (Kosswig 1955), where there are no effective barriers to prevent gene flow between populations.TR3 corresponded to the O * haplotype (see Table 2 and Frati et al. 1998), which is distributed in Bulgaria.Furthermore, TR7 corresponded to the C13 haplotype (see Table 2 and Inoue et al. 2007), which is distributed in Hokkaido, Japan.
Based on three different data sets of the mitochondrial Cyt b sequences (250 bp, 354 bp and 375 bp), network analyses suggested that the majority of Turkish red fox, including 49 samples, belonged to Eurasia phylogroup, and that two haplotypes (TR7 and TR10) from SW Anatolia, Turkey clustered together with three haplotypes from Hokkaido, Japan (Inoue et al. 2007) (Figs. 2, 3).The existence of different haplotypes within Turkey may be due to the maternal gene flow, a founder effect or mutations accumulated throughout time.The Eurasian haplotypes in the Asian part of Turkey (Anatolia) are a signal of maternal gene flow via the European part of Turkey (Turkish Thrace) in the past.As stated by many authors (taberlet et al. 1998, Hewitt 1999, Inoue et al. 2007, Aubry et al. 2009), this admixture of Eurasian haplotypes and local haplotypes in Anatolia may be a result from dispersal or migration of red foxes during the glacial and interglacial periods of the Pleistocene.However, the close relationship between red fox haplotypes from   with previous studies (Inoue et al. 2007, Aubry et al. 2009, Oishi et al. 2011, Yu et al. 2012a).The occurrence of two distinct phylogroups was reported by Inoue et al. (2007) and Oishi et al. (2011)  There are many studies that provided useful data which assists to understand the phylogeographic patterns of numerous species across Turkey and the surrounding regions of the Mediterranean Basin.These publications cover various species like red deer (Ludt et al. 2004), brown trout (Bardakci et al. 2006), reptiles (Joger et al. 2007), freeliving mice (Macholan et al. 2007), house mouse (Rajabi-Maham et al. 2008), ground squirrels (Gündüz et al. 2007, Kryštufek et al. 2009), Middle Eastern tree frogs (Gvozdik et al. 2010), voles (Thanou et al. 2012), orthoptera (Çiplak et al. 2005, Kaya et al. 2013) and herptiles (Poulakakis et al. 2013) from Holoarctic region that used to study the historical processes influencing their current distribution.Additionally, it has been suggested that the modern red fox (V.vulpes) originated from Asia Minor (Anatolia) or North Africa around the Mediterranean Basin (Hersteinsson & Macdonald 1982, Szuma 2008), or China (Kurtén 1968).However, the origin of the modern red fox is discussed controversially, since the North African and Turkish red foxes are not investigated in detail thus far.Although we aimed to add some insights into the current discussion, our study has not showed a sufficient evidence that the modern red fox originated from Asia Minor (Anatolia) due to the clustering positions of Turkish haplotypes in all analyses (Figs. 2 -5).
When it comes to the remarkable clustering of two haplotypes from SW Anatolia with haplotypes from Hokkaido, the indicated close relationship between haplotypes from the both areas is hardly comprehensible from the zoogeographic point of view, and it may be the result of independent evolution (random homoplasy).Given that the level of haplotype diversity in Turkey and the position of TR9 in haplotype networks (Figs. 2, 3), the northern part of Turkey, where the Turkish haplotype TR9 is distributed, may be considered as the Pleistocene refugium for red fox (Table 1, Fig. 1).The north Anatolian part of Turkey was also accepted to be a possible refugium for several species (Santucci et al. 1998, Hewitt 1999, 2000, Wielstra et al. 2013).In contrast, recent studies, not including Turkey, suggest that red foxes in Eurasia survived the last glacial maximum as a single, large, interbreeding population (Teacher et al. 2011, Edwards et al. 2012, Kutschera et al. 2013).
In conclusion, compared with other parts of the red fox's range, the occurrence of at least two main phylogenetic groups and more than one subphylogroups (Figs. 2 -5) and the finding of high genetic variation in red foxes from Anatolia confirmed the importance of the studied region.In this sense, to elucidate in depth the Pleistocene zoogeography of the Turkish red fox and the phylogeographic pattern of the red fox through-out the species range, one should use much larger data sets from mitochondrial DNA and nuclear DNA markers and, in the meantime, collect more red fox samples from the other locations of Turkey, North Africa, North-east Europe, Caucasus and the Central Asia as well as South Asia including Iran, Southern China, Pakistan, northern India.Additionally, samples from Australia, introduced in the 1800's, may also add data sets since mitochondrial D-loop sequences are deposited in GenBank.

Fig. 2 .
Fig. 2. Median-joining network reconstructed from haplotypes of Vulpes vulpes in the North Hemisphere, based on 250 bp sequences of the mitochondrial cytochrome b.Bold indicates the Turkish haplotypes.See Tables 1, 2 for the sample abbreviations.

Fig. 3 .
Fig. 3. Median-joining network reconstructed from haplotypes of Vulpes vulpes in the North Hemisphere, based on 375 bp sequences of the mitochondrial cytochrome b.Bold indicates the Turkish haplotypes.See Tables 1, 2 for the sample abbreviations.

Fig. 5 .
Fig. 5. Bayesian phylogenetic tree reconstructed from haplotypes of Vulpes vulpes in Eurasia, based on 375 bp sequences of the mitochondrial cytochrome b.Numbers above branches show the Bayesian posterior probabilities.Bold indicates the Turkish haplotypes.See Tables 1, 2 for the sample abbreviations.

Fig. 4 .
Fig. 4. Bayesian phylogenetic tree reconstructed from haplotypes of Vulpes vulpes in the North Hemisphere, based on 250 bp sequences of the mitochondrial cytochrome b.Numbers above branches show the Bayesian posterior probabilities.Bold indicates the Turkish haplotypes.See Tables 1, 2 for the sample abbreviations.
using mitochondrial DNA (Cyt b and control region) and microsatellites in Japanese red fox, respectively, by Aubry et al. (2009) using mitochondrial DNA (Cyt b and control region), by Yu et al. (2012a) using the cyt b in Korean red fox.

Table 1 .
). List of Turkish samples and sequences/haplotypes used in the present study (* Map numbers are locality codes in Figure1).

Table 2 .
The sample abbreviations for haplotype information of the previous reports used in median-joining networks and Bayesian trees.