Email updates

Keep up to date with the latest news and content from BMC Evolutionary Biology and BioMed Central.

Open Access Highly Accessed Research article

A range-wide synthesis and timeline for phylogeographic events in the red fox (Vulpes vulpes)

Verena E Kutschera1*, Nicolas Lecomte23, Axel Janke14, Nuria Selva5, Alexander A Sokolov6, Timm Haun1, Katharina Steyer7, Carsten Nowak7 and Frank Hailer1*

Author Affiliations

1 Biodiversity and Climate Research Centre (BiK-F), Senckenberg Gesellschaft für Naturforschung, Senckenberganlage 25, 60325 Frankfurt am Main, Germany

2 Department of Arctic and Marine Biology, University of Tromsø, N-9037 Tromsø, Norway

3 Department of Environment, Government of Nunavut, X0A0L0 Igloolik, Canada

4 Institute for Ecology, Evolution and Diversity, Goethe University Frankfurt, Max-von-Laue-Straße 13, 60438 Frankfurt am Main, Germany

5 Institute of Nature Conservation, Polish Academy of Sciences, Mickiewicza 33, 31-120 Krakow, Poland

6 Ecological Research Station of the Institute of Plant and Animal Ecology, Russian Academy of Sciences, Labytnangi, 629400, Russia

7 Conservation Genetics Group, Senckenberg Research Institute and Natural History Museum Frankfurt, Clamecystraße 12, 63571 Gelnhausen, Germany

For all author emails, please log on.

BMC Evolutionary Biology 2013, 13:114  doi:10.1186/1471-2148-13-114

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2148/13/114


Received:1 February 2013
Accepted:29 May 2013
Published:5 June 2013

© 2013 Kutschera et al.; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Many boreo-temperate mammals have a Pleistocene fossil record throughout Eurasia and North America, but only few have a contemporary distribution that spans this large area. Examples of Holarctic-distributed carnivores are the brown bear, grey wolf, and red fox, all three ecological generalists with large dispersal capacity and a high adaptive flexibility. While the two former have been examined extensively across their ranges, no phylogeographic study of the red fox has been conducted across its entire Holarctic range. Moreover, no study included samples from central Asia, leaving a large sampling gap in the middle of the Eurasian landmass.

Results

Here we provide the first mitochondrial DNA sequence data of red foxes from central Asia (Siberia), and new sequences from several European populations. In a range-wide synthesis of 729 red fox mitochondrial control region sequences, including 677 previously published and 52 newly obtained sequences, this manuscript describes the pattern and timing of major phylogeographic events in red foxes, using a Bayesian coalescence approach with multiple fossil tip and root calibration points. In a 335 bp alignment we found in total 175 unique haplotypes. All newly sequenced individuals belonged to the previously described Holarctic lineage. Our analyses confirmed the presence of three Nearctic- and two Japan-restricted lineages that were formed since the Mid/Late Pleistocene.

Conclusions

The phylogeographic history of red foxes is highly similar to that previously described for grey wolves and brown bears, indicating that climatic fluctuations and habitat changes since the Pleistocene had similar effects on these highly mobile generalist species. All three species originally diversified in Eurasia and later colonized North America and Japan. North American lineages persisted through the last glacial maximum south of the ice sheets, meeting more recent colonizers from Beringia during postglacial expansion into the northern Nearctic. Both brown bears and red foxes colonized Japan’s northern island Hokkaido at least three times, all lineages being most closely related to different mainland lineages. Red foxes, grey wolves, and brown bears thus represent an interesting case where species that occupy similar ecological niches also exhibit similar phylogeographic histories.

Keywords:
Carnivores; Divergence time estimate; Generalist; mtDNA control region; Phylogeography; Vulpes

Background

While current population genetic structuring tends to be weak in arctic mammalian specialists [1,2], species in temperate regions usually show more pronounced structuring [3-5], due to their survival in different refugia during the last glacial maximum (LGM) (reviewed in [6,7]). An exception to this trend is observed in some temperate-zone species with generalist habitat requirements and high mobility, like large carnivores that occur across a wide habitat and climatic gradient. Indeed, grey wolves (Canis lupus) and brown bears (Ursus arctos) show an overall weak phylogeographic structure with several widely distributed lineages [8-15].

Similar to the brown bear and the grey wolf, the red fox (Vulpes vulpes) is distributed across all northern continents (Europe, Asia, and North America), being the most widely distributed carnivore in the world [16]. The high mobility and adaptability of the red fox to different habitats and climates is reflected by its earliest appearance in the fossil record outside the southern refugia shortly after the last glacial maximum (LGM). Some 13,500 years before present, for instance, the red fox re-appeared in Northern Germany close to the ice sheets [17]. Daily distances of more than 10 km are common [18,19], and the longest recorded distance covered by a red fox was 302 km within less than a year’s time [20].

Although red fox phylogeography using DNA sequence data from a ca. 268 – 342 bp fragment of the mitochondrial (mt) control region has been investigated before [21-27], these studies had a regional focus or sampling gaps in Siberia and Asia (Figure  1). In addition, so far the published data has not been collated to investigate range-wide processes, and no timeline has yet been established for phylogeographic events. Previous studies described two major red fox lineages – one with a Holarctic distribution and a Nearctic lineage consisting of three sublineages (widespread lineage, eastern lineage, mountain lineage) [23,24,27]. A study on Japanese red foxes found two main lineages in Japan: one lineage that was exclusively found on Japan’s northern island Hokkaido (Hokkaido II), and another lineage that comprised three Japanese sublineages occurring on Hokkaido and on Japan’s main southern islands Honshu and Kyushu (Hokkaido Ia, Hokkaido Ib, Honshu/Kyushu) along with mainland Asian red foxes [21].

thumbnailFigure 1. Map of sample locations for all sequences used in this study (published data and newly generated data). Current distribution range of the red fox [84] is shown in light grey. Black stars with white numbers indicate sampling locations for new data generated in this study. Circles indicate sampling regions from previous studies, with black numbers denoting numbers of published sequences for these regions. Details on all sequences used in this study are given in Additional file 1.

Additional file 1. Sequence information for all individuals analyzed in this study. Individual IDs, haplotype frequencies, consecutive haplotype numbers, individuals with identical sequence, geographic origins, age assumed for the BEAST analyses, GenBank accession numbers, corresponding abbreviations used in the original source studies, and the corresponding references are provided. The designation of individuals is explained in detail in the table caption.

Format: PDF Size: 339KB Download file

This file can be viewed with: Adobe Acrobat ReaderOpen Data

We here present novel data from red foxes of various Eurasian populations, notably including individuals from central Siberia, a biogeographically important region with the potential to link European with East Asian and/or North American populations. In a range-wide synthesis of published and publicly available control region sequences combined with newly generated data (Figure  1), the identity and geographic distribution of previously described lineages is validated. Based on this novel assessment of mtDNA structure in red foxes, we establish a timeline of major phylogeographic events using a Bayesian coalescence approach with multiple fossil tip and root calibration points. We compare these results for red foxes to previously published findings from other carnivores with a Holarctic distribution, allowing us to identify common phylogeographic patterns and processes.

Results

In our analysis of 729 red foxes (a map with the geographic locations is given in Figure  1 and a complete compilation is found in Additional file 1), we found 95 variable sites that defined 175 haplotypes in a 335 bp-alignment of the mitochondrial control region (Table  1). Haplotype diversity for the whole dataset was 0.948 +/− 0.005, and nucleotide diversity was 0.057 +/− 0.028 (Table  1). The 52 newly sequenced individuals from Siberia, Germany, Poland, and Finland formed 25 haplotypes. Most of these haplotypes (n = 22) had not been encountered in previous studies (Additional file 2). All newly obtained sequences have been submitted to the EMBL database [EMBL:HF677203-HF677255].

Table 1. Summary statistics of genetic variability of major red fox mtDNA control region lineages

Additional file 2. Sample size and number of mtDNA control region haplotypes for newly sequenced red foxes. This pdf-file contains a table giving details on haplotypes (novelty, haplotype-sharing) that were reconstructed from 52 newly sequenced red foxes.

Format: PDF Size: 50KB Download file

This file can be viewed with: Adobe Acrobat ReaderOpen Data

Major mitochondrial clades and population structuring

We confirmed all previously described lineages (Figures  2 and 3, Additional file 3). All Nearctic lineages and the Japanese Hokkaido II and Honshu/Kyushu lineages were distinct in a median-joining network (Figure  3), although some received less than 95% posterior support in the BEAST analysis (Figure  2). However, the most basal split within each of these lineages received high posterior support in the BEAST analysis (Figure  2). A phylogenetic analysis of haplotype data conducted in MrBayes recovered a tree with a topology congruent with the tree obtained from BEAST (Additional file 3), in accordance with previous lineage definitions. Therefore, despite some uncertainty with regard to their placement in the red fox mtDNA phylogeny, those regional lineages represent distinct clades that capture aspects of the evolutionary history of red foxes. Our analysis and discussion of regional lineages therefore focuses on groupings with high statistical support and/or lineages that were specifically defined in previous studies.

thumbnailFigure 2. Maximum clade credibility tree with time estimates for colonization events and basal diversifications within red fox mtDNA control region lineages. For the dating of phylogeographic events we used a combined approach, utilizing tip dates based on published ancient DNA red fox sequences [25], plus the arctic fox as exterior calibration point. This tree shows the results for a root height prior of 1.75-4 million years (My), according to the 95% credibility interval in Perini et al. [80] for the divergence time between red and arctic fox (Table  2). White: Nearctic lineage haplotypes; grey: Japanese haplotypes (lineages Honshu/Kyushu and Hokkaido II); black: Holarctic lineage haplotypes, including Japanese lineages Hokkaido Ia and Ib; ka: thousand years. Nodes marked with an asterisk were supported by posterior probability values >0.95. Samples used for tip calibration are marked with a ♦ symbol. Median ages and 95% highest posterior density ranges in brackets show the estimated ages of major lineages, and of the most basal nodes within these lineages. Our discussion focuses on lineages/nodes with ≥0.95 statistical support, recognizing that longer mtDNA sequences will be required to resolve larger proportions of the red fox mitochondrial phylogeny (see [58-62]).

thumbnailFigure 3. Median joining network of genetic variation at a 335 bp fragment of the mitochondrial control region in 693 red foxes. 35 ancient samples from Teacher et al. [25] and one modern sample from Valière et al. [28] were excluded due to their shorter sequence length. White: Nearctic lineage haplotypes; grey: Japanese haplotypes (all Japanese lineages); black: Holarctic lineage haplotypes, except for North American, Japanese and central Siberian haplotypes; black with white circle: North American haplotypes (Holarctic lineage); black with white stripes: central Siberian haplotypes (Holarctic lineage).

Additional file 3. Bayesian inference tree of red fox mtDNA control region haplotypes. This png-file contains a Bayesian inference tree that was based on 175 haplotypes reconstructed in MrBayes. All major lineages are indicated by square brackets. Interesting haplotypes within the Holarctic lineage are indicated as follows: grey stars: Japanese Hokkaido Ia and Ib haplotypes (Holarctic lineage); black stars: Serbian haplotypes; black stars with parallel white stripes: central Siberian haplotypes; black stars with white edge: North American haplotypes (Holarctic lineage).

Format: PNG Size: 1.5MB Download fileOpen Data

Nucleotide and haplotype diversity were highest for geographically widely distributed lineages (Holarctic, Nearctic, and Nearctic widespread lineages), and lowest for regionally restricted lineages (Japanese Hokkaido II and Honshu/Kyushu lineages, Nearctic mountain and eastern lineages) (Table  1).

North American red foxes fell into four mitochondrial groups: a Nearctic lineage that comprised three sublineages (mountain, eastern, and widespread lineages), and a Holarctic lineage that was composed of North American and Eurasian red foxes (Figures  2 and 3). These lineages were previously described and discussed in detail by Aubry et al. [23], Sacks et al. [24], and Statham et al. [27].

All Eurasian red foxes were placed in the Holarctic lineage (Figures  2 and 3). In contrast to the Nearctic lineage, most support values for phylogenetic groupings within the Holarctic lineage were low: no geographically restricted sublineages in the Holarctic lineage received high statistical support, except for two sublineages from Japan (Figures  2 and 3; and see below). Haplotypes from other geographic regions did not form monophyletic groups, but were scattered across the Holarctic lineage (Figures  2 and 3). However, the few instances of haplotype sharing occurred only between geographically neighboring populations (e.g., Germany, Switzerland, France; see Additional file 1). This weak phylogeographic structuring within the Holarctic lineage was consistent with previous findings by Teacher et al. [25] and Edwards et al. [26], who discussed red fox phylogeographic patterns in Europe more in detail.

All 52 newly sequenced red foxes from Siberia, Germany, Finland, and Poland belonged to the Holarctic lineage (Figures  2 and 3). The six central Siberian individuals formed three not previously described haplotypes that were not especially closely related to each other (Figure  3, Additional files 1 and 2). Also the Finnish red fox carried a novel haplotype (F01). The thirty newly sequenced German individuals formed ten haplotypes. Seven of these were new; three haplotypes were identical to previously published sequences from French, Swiss or other European populations [23,25,28] (Additional files 1 and 2). One German haplotype (individuals D408, D655, and D660) was shared with a Polish red fox (POLI68) from the present study. The newly sequenced Polish red foxes formed thirteen haplotypes (Figure  3, Additional files 1 and 2), twelve of which had not been described previously. Besides the haplotype described above (POLI68), another was identical to a German haplotype from the present study (D08) and to a Swiss haplotype from a previous study (Swit12) [25,28] (Additional file 1). One Polish haplotype was closely related to a Serbian haplotype (I) that was previously described as being distinct from all other Serbian red foxes [22] (Serbia 3; see Figure  3).

The other Serbian red foxes clustered separately into two groups that were distinct from each other (Serbia 1 and Serbia 2) within the Holarctic lineage, confirming Kirschning et al. [22]. In contrast to Kirschning et al. [22] findings, however, the two groups were separated by several Eurasian haplotypes and distinct from the rest of the Holarctic diversity in the median-joining network (Figure  3). In the Bayesian inference tree (Additional file 3), all three Serbian lineages received less than 95% posterior support.

The Japanese samples formed four distinct groups that fell within the Holarctic diversity (Figures  2 and 3, Additional file 3). These Japanese lineages were separated by intermediate haplotypes found on the Eurasian mainland (Figure  3). Unlike in previous analyses [21,29], our extended sampling and range-wide synthesis revealed that the Honshu/Kyushu lineage was not closely related to Hokkaido Ia, Hokkaido Ib, and Asian mainland red foxes, but formed a distinct lineage (Figures  2 and 3, Additional file 3). Red foxes from the northern island Hokkaido formed three separate groups. The Hokkaido II lineage remained clearly distinct from the other red foxes, as originally described by Inoue et al. [21] and confirmed using mitochondrial cytochrome b data by Yu et al. [29] (Figures  2 and 3). Hokkaido I red foxes clustered into two subgroups in the Holarctic lineage (Figures  2 and 3) consistent with Inoue et al. [21] findings of sublineages Hokkaido Ia and Hokkaido Ib. In contrast, however, our range-wide synthesis recovered Eurasian mainland haplotypes that separated the Hokkaido Ia and Ib groups from each other (Figure  3). Further, one haplotype from Hokkaido Ib (JH9) was identical to a haplotype found by Aubry et al. [23] on the Asian mainland (AS70, see Figure  3 and Additional file 1), supporting the close relatedness between these populations.

Timing of phylogeographic events in the red fox

The influence of different root heights on the divergence time estimate for the red fox from the arctic fox (Vulpes lagopus) was evaluated by performing different simulations in BEAST. However, each BEAST run converged on the youngest time frame allowed by the prior for the speciation event (i.e., close to the different minimum root heights of 1.75 million years (My) and 5.1 My, respectively; see Table  2). Similarly, a BEAST run that only applied tip dating (i.e., with no constraint on minimum root height) yielded a divergence time estimate for red/arctic foxes of 301 (101–611) thousand years (ka). This time estimate is less than the first appearance of the red fox in the fossil record, ca. 0.5 – 1 million years ago (Mya) [30,31], reinforcing the value of mixed tip/root calibration approaches in BEAST [13,32]. Hence, our data do not allow an accurate assessment of the divergence time between red and arctic foxes - the focus of this paper is on much more recent phylogeographic events within red foxes.

Table 2. Comparison of BEAST dating results employing a relaxed-clock approach with combined tip (interior) and root (exterior) calibration

Importantly, our analyses in BEAST indicated that this uncertainty about the root height of the tree (i.e., the divergence from arctic foxes) did not have a major effect on the dating of evolutionary events within red foxes. Comparing the BEAST simulations with a root height of minimum 0.5 My to the simulation setting it to at least 5.1 My – the two most extreme scenarios including a root height and tip dating – the estimated divergence time for the red/arctic fox speciation event varied by a factor of 10, whereas the median time estimates for phylogeographic events within red foxes varied only by a factor of 1.4-3, with overlapping confidence intervals (Table  2). It thus appeared that the inconsistent root height (red/arctic fox divergence time) only slightly impacted our time estimates for phylogeographic events within red foxes, which is likely in part due to our additional use of interior tip calibrations based on ancient DNA sequences from known-age fossil remains.

BEAST runs yielded posterior substitution rate estimates of 33.2% for Scenario 1 (with a root height of 1.75 My), and 25.6% or 41.9% for Scenarios 2 and 3, respectively (Table  2). A recent study of red fox mtDNA obtained a similar mutation rate estimate (ca. 26.29-33.81%, depending on whether the substitution model included a gamma correction; Edwards et al. [26]). It is interesting to note that these results [26] were based on a partial control region fragment combined with a partial cytochrome b fragment, precluding direct comparisons of the estimates. Similar reasoning applies to the mutation rate estimate (28.8%) obtained for mtDNA in brown and cave bears [33], that is based on a fragment that only overlaps partially with our alignment.

Due to low posterior support values for many internal nodes (Figure  2 and Additional file 3), our analyses do not allow to identify the next most closely related sequence of some of the described lineages. However, regardless of their exact phylogenetic placement, estimated divergence times of regional lineages from their most closely related sequence remained largely constant (Table  2).

Based on the most basal divergences among extant lineages (Figure  2, Table  2), red foxes started to diversify in Eurasia during the end of the Mid Pleistocene (Figure  4). North America was colonized independently by several lineages from this Holarctic diversity (Figures  2 and 4, Table  2), as indicated by simulations in BEAST. The colonization event forming the Nearctic lineage happened around the Mid or Late Pleistocene (Figures  2 and 4, Table  2). All three Nearctic sublineages (eastern, mountain, and widespread) were formed before the LGM (Figures  2 and 4, Table  2). The close relationship between North American and Eurasian Holarctic lineage red foxes (Figures  2 and 3) indicated that the North American Holarctic lineage colonized North America much more recently than the Nearctic lineage, probably around the LGM (Figure  4).

thumbnailFigure 4. Map outlining major phylogeographic events in red foxes as reconstructed using a Bayesian coalescence approach with multiple fossil tip and root calibration points. Current distribution range [84] is shown in foxy red. White: Nearctic lineage; grey: Japanese lineages; black: Holarctic lineage (excluding Japan-restricted lineages). MP: Mid Pleistocene, LP: Late Pleistocene, LGM: Last Glacial Maximum, Hol: Holocene. Note that not all sublineages within the Holarctic lineage are currently distributed across Eurasia and North America; only some lineages show extensive range expansions (Figures  2 and 3, and Additional file 3).

Also Japan was colonized several times independently by individuals from the Holarctic lineage. The southern islands Honshu and Kyushu were colonized during the Late Pleistocene, or around the LGM (Figures  2 and 4, Table  2). The northern island Hokkaido was colonized several times. The lineage Hokkaido II arrived most likely during the end of the Late Pleistocene (Figures  2 and 4, Table  2). Due to the close relationship of Hokkaido Ia and Hokkaido Ib red foxes to Eurasian mainland red foxes (Figures  2 and 3), the two groups most likely colonized Hokkaido more recently than the Hokkaido II lineage did, probably around the LGM (Figure  4).

Demography

When we tested the exponential population growth model in BEAST, the exponential growth rate fluctuated around zero, so we could not reject the constant population size model for the entire dataset. Fu’s FS showed a signal of population growth across all red fox sequences, with a highly significant (p ≤ 0.001) value of −23.363 (Table  1). Except for the Nearctic widespread lineage, all major lineages had a significant (p ≤ 0.02) FS, with negative values ranging from −0.189 for the Japanese Hokkaido II lineage to −23.701 for the Holarctic lineage (Table  1). Star-like structures in a median joining network (Figure  3) also indicated evolutionarily recent population growth for these lineages.

We calculated mismatch distributions for all lineages in Arlequin. For the Nearctic lineage, the Nearctic eastern and widespread lineages, and the Japanese Hokkaido II lineage, population growth was confirmed by mismatch distribution analyses, where simulations in Arlequin did not differ significantly (p > 0.05) from expectations under the sudden expansion model (Table  3). The Nearctic lineage started to diversify around the Late Pleistocene/LGM, as further supported by our dating in BEAST of the most basal bifurcation within each lineage (Figure  2, Table 2 and  3). The estimated sudden expansion for the Nearctic eastern lineage of 23 (11–35) ka (assuming 7.1% per-lineage substitution rate per My), and 5 (2–7) ka (assuming 33.2% per-lineage substitution rate per My) (Table  3) overlapped with the period indicated by the BEAST analyses (Figure  2, Table  2). Population expansion of the Japanese Hokkaido II lineage started about 36 (0–85) ka with 7.1% per-lineage substitution rate per My, and about 8 (0–18) ka with 33.2% per-lineage substitution rate per My (Table  3), similar to the date estimate from the BEAST analysis of 6 (1–15) ka (Figure  2). The very recent date estimates for the most basal bifurcation in Hokkaido II during the Holocene (Figure  2, Table  3) were possibly a result of genetic drift removing more ancient lineages in this clade (e.g., [34]). Because deviation from the sudden expansion model in Arlequin was significant for the remaining major clades, we could not determine the timing for any onset of population expansion for them. However, based on the BEAST results, the Honshu/Kyushu population diversity was slightly older than the Hokkaido II population diversity, with 95% credibility intervals spanning the LGM (Figure  2, Table  3).

Table 3. Mismatch distribution analysis under a sudden expansion model and time since expansion calculated for different mitochondrial lineages

The two methods we used to detect population growth are characterized by different levels of sensitivity: tests based on mismatch distribution analyses are conservative and less powerful than Fu’s FS to detect population growth [35]. This might explain inconsistencies in our demographic inferences: Fu’s FS was negative and significant for the Nearctic mountain lineage and the Japanese Honshu/Kyushu lineage (Table  1), whereas mismatch distribution analyses failed to detect population growth (thus not enabling us to obtain an estimate of the timing from the mismatch distribution analysis for these lineages). The Nearctic widespread lineage showed the opposite result: while mismatch distribution analyses indicated population growth (Table  3), Fu’s FS was positive (Table  1). Our inferred timing of the onset of sudden expansion of this lineage from mismatch distribution analysis (about 280 (15–522) ka with 7.1% per-lineage substitution rate per My, and about 60 (3–112) ka with 33.2% per-lineage substitution rate per My; Table  3) should therefore be taken with caution, especially since the time estimate inferred from BEAST (24 (9–46) ka; Figure  2; see also Table  2) was younger.

Discussion

Our synthesis of 677 published and publicly available sequences together with 52 newly obtained sequences includes previously unsampled geographic regions (e.g., central Siberia; Figure  1), and confirms previous classifications of mitochondrial lineages in red foxes [21,23,24,27]. In this study, we delineate the range-wide timing and pattern of phylogeographic events for this widespread carnivore. During the Mid Pleistocene, a period characterized by repeated climatic oscillations, red foxes started to diversify. One lineage (the Holarctic lineage) today occurs across most of the entire distribution range, including North America, Europe, and Asia. In contrast, other red fox lineages are regionally restricted. During the Late Pleistocene and Holocene, North America and Japan were colonized several times independently by red foxes (Figures  2 and 4), likely at times when landbridges connected these islands to adjacent landmasses (except for Honshu/Kyushu, see below). Our dating of subsequent diversification events suggests that demographic expansions in many red fox populations occurred since the Late Pleistocene. Our suggested timeline and pattern of phylogeographic events in red foxes closely resembles the scenarios described for the ecologically similar and co-distributed grey wolves and brown bears [8,10,13,36], reinforcing that ecological and climatic factors had similar effects on temperate zone species.

Phylogeographic history of red foxes

Mid Pleistocene

According to the fossil record [37], red foxes were already present in Eurasia during the Mid Pleistocene, around 300 ka, a timing consistent with our finding of the most basal diversification within red foxes during that time (Figure  2). For extended time periods since the Mid Pleistocene, North America was connected to Eurasia via the Bering landbridge [38], at glacial times of lower sea level, allowing species like the red fox to colonize North America several times independently. Red fox fossil remains from Alaska indicate that North America was first colonized during the Illinoian glaciation (ca. 300–130 ka) [39]. This period overlaps with our estimate for the emergence of the Nearctic lineage (Scenario 2, Table  2), likely corresponding to the colonization of North America from Eurasia around the end of the Mid Pleistocene (Figures  2 and 4).

Late Pleistocene

The fossil record indicates that some North American red foxes persisted through the Late Pleistocene glaciations (Wisconsin; ca. 100–10 ka) south of the ice sheets [40-42]. It would be interesting to use ancient DNA techniques to verify whether such remains indeed belong to the Nearctic red fox mitochondrial lineage. Consistent across a broad range of different rate calibrations (Table  2), we found a Late Pleistocene diversification of the Nearctic lineage into three sublineages (Figures  2 and 4), the eastern, mountain, and the widespread lineages.

During the same time frame, the Holarctic lineage further diversified in Eurasia (Figures  2 and 4). From this diversity, Japan was colonized several times independently. Japan’s main southern islands Honshu and Kyushu have been isolated from the Eurasian mainland and Hokkaido in the north since the Mid Pleistocene (reviewed in [43]). The Tsugaru Strait separating Hokkaido from Honshu/Kyushu represents a biogeographic barrier (Blakiston’s line [44]) for many species since the Mid Pleistocene (reviewed in [43]), including the red fox. However, red foxes appear to have colonized Honshu and Kyushu during the Late Pleistocene (Figures  2 and 4). The European fossil record indicates potential human influence on red fox dispersal (based on the common finding of red fox remains in European archaeological assemblages; [17]). Humans first reached Japan ca. 50 ka (reviewed in [45]), consistent with our red fox colonization time estimate for Honshu and Kyushu. Red foxes might thus have reached these islands through human introduction, although non-human facilitated rafting or temporary ice bridges connecting Honshu and Hokkaido [46] are also plausible alternatives.

During most of the Late Pleistocene glacial phases, Hokkaido was repeatedly connected to the Eurasian mainland via the Russian island Sakhalin [43]. Indeed, red foxes from the Hokkaido II lineage colonized Hokkaido during the Late Pleistocene (Figures  2 and 4), most likely from the mainland via these northern landbridges.

Around the last glacial maximum (LGM)

The phylogeographic structure found in North American red foxes has been interpreted to result from range fragmentation by ice sheets during glacial maxima [23]. In contrast, we did not find a strong signal of red fox survival during the LGM in distinct and isolated southern refugia in Eurasia. In fact, the Holarctic lineage appears to be geographically largely unstructured and shows only weak signs of past population fragmentation (Figures  2 and 3; see also [25,26]). This could support the notion of red foxes as habitat generalists that were able to survive the LGM period in the vast arctic steppes, as shown for a now-extinct wolf ecotype [34]. However, red fox remains have not been found in the arctic steppe regions of northern and central Europe during the LGM [17]. Indeed, contemporary red fox distribution and population ecological studies suggest that the species has a northern (and upper latitudinal) limit in its distribution, where more arctic-adapted species like the arctic fox have a competitive advantage (e.g., [47]). Survival of red foxes during the LGM in southern refugia is supported by several lines of evidence: (i) our finding of (albeit weakly) differentiated mitochondrial sublineages in Serbia (Figure  3), (ii) genetic differentiation of Iberian red foxes from other European populations [26], (iii) signs of low connectivity among different Mediterranean populations [48], and (iv) the restriction of red foxes to southern European refugial regions during the LGM [17]. Similar to findings from brown bears [10], the high dispersal capability of red foxes and their rapid recolonization of suitable habitats after deglaciation may have led to gene flow among refugia, preventing the development of pronounced phylogeographic structuring. Our finding of signals of population expansions in many red fox lineages (Tables  1 and 3) likely reflects this postglacial demographic growth.

The following conditions may therefore have prevented the formation of deeply separated mitochondrial lineages in Eurasian red foxes: only short periods (ca. 10 ka [17]) of geographic restriction in refugia for genetic differentiation, and potential for some gene flow among such refugial regions in periods of temporary warming within longer climatic cold phases. Ongoing range expansion of red foxes north and into higher altitudes into traditional arctic fox habitats (e.g., [49-51]) may mirror this situation of rapid recolonization of northern habitats.

The extensive mitochondrial gene flow among Eurasian red fox populations also reached parts of North America, forming a vast Holarctic population (Figure  4). Ice-free regions of Alaska and the Yukon (Beringia) were connected to Eurasia via the Bering landbridge during glacial maxima, but separated from regions south of the Laurentide and Cordilleran ice sheets [7]. Due to postglacial sea level rise, the Holarctic North American population was eventually isolated from the rest of the Holarctic lineage when the Bering landbridge was closed. Today, North American and Eurasian haplotypes from the Holarctic lineage are still intermingled (Figures  2 and 3), confirming their recent evolutionary separation. However, we found no haplotype sharing between North American and Eurasian red foxes (but note the shared haplotype between mainland Asia and Hokkaido; see Additional file 1), likely reflecting a post-LGM interruption of trans-Beringian gene flow. Similarly, the British Isles belonged to the Eurasian landmass during much of the Late Pleistocene and early Holocene. Edwards et al. [26] found red foxes from the British Isles to be only weakly differentiated from European main land foxes, resulting from recent postglacial isolation. Japan’s northern island Hokkaido was also repeatedly part of the Eurasian landmass during the Late Pleistocene, which is reflected by the presence of at least three distinct red fox lineages (Hokkaido Ia, Hokkaido Ib, and Hokkaido II; Figures  2 and 3; [21]). Hokkaido I haplotypes (Ia and Ib) are closely related to Eurasian mainland haplotypes (Figures  2 and 3), suggesting that they were isolated only recently from the rest of the Holarctic lineage. This may have occurred after the LGM (Figures  2 and 4) when rising sea levels isolated Hokkaido from the Asian mainland [43].

Genetic and paleontological data indicate that population expansion after the LGM occurred rapidly in European red foxes ( [17,26]; see also the signals of demographic expansion in Table  3). The finding of sublineages within the Holarctic diversity in southern Europe [26,48] indicates that some southern red fox populations contributed less to the postglacial recolonization of Eurasia, while other lineages showed wide-ranging dispersal, even across Beringia (Figure  4).

Despite these signals of extensive wide-ranging gene flow in red foxes, mtDNA also shows some signals of currently restricted gene flow, even within continents. Studies reporting local adaptations by red foxes in some regions that are discussed below confirm this view. In our dataset, the only instances of haplotype sharing occur among neighboring populations in Europe, and a single case involving Hokkaido and the Asian mainland (Additional file 1). Similarly, studies on red foxes using nuclear microsatellites have shown dispersal restrictions on shorter time scales [24,52-54], consistent with lower levels of gene flow today.

Common phylogeographic trends in Holarctic carnivores

Two other large carnivores that have a Holarctic distribution and generalist habitat requirements are the brown bear and the grey wolf. Both species are highly mobile and flexible regarding their habitat requirements. Besides their ecology, they share remarkably similar phylogeographic patterns with the red fox. All three species are characterized by very widespread Holarctic mitochondrial lineages that are distributed across Eurasia and North America, and only some locally restricted lineages - indicating dispersal limitations in some regions of the Holarctic [8,9,11-13,34,55].

Independent colonization events from Eurasia led to the establishment of several endemic North American lineages such as grey wolves [34] and brown bears [13]. Those lineages survived the Late Pleistocene glaciations south of the North American ice sheets. Since the LGM period, additional Eurasian lineages immigrated via the Bering landbridge into the Nearctic, leading to current denomination of those lineages as Holarctic (Figure  4). Brown bears also inhabit the Japanese island Hokkaido. Similar to the situation in red foxes, Hokkaido was colonized at least three times by brown bears. Hokkaido harbors three distinct brown bear lineages, each being most closely related to different mainland lineages [56].

The grey wolf is another extensively studied example of a Holarctic generalist that displays a worldwide weak phylogeographic pattern, but with some geographically restricted mitochondrial lineages. It appears that especially some southern wolf populations did not expand after the Pleistocene glaciations: two distinct lineages survived south of the Himalayas [9]. In brown bears, a similar situation has been described for Syria and Iran [11,13], and in red foxes analogous evidence exists for Iberia [26,48] and Serbia (Figure  3; data from Kirschning et al. [22]).

Likely reflecting even more recent processes, grey wolves have been shown to be locally adapted to specific habitat and foraging conditions in Pacific temperate coastal rainforests [12,15]. Further, eastern European grey wolf population structure appears to be correlated with ecological factors [36,57]. As reviewed by Sacks et al. [24], North American montane red fox populations show physiological and morphological adaptations to cold climate, and are genetically distinct from other red fox populations (Nearctic mountain lineage) [23]. A refined sampling in this region supported an indigenous origin of the Sacramento Valley population, which differs in body size from the montane red fox populations [24].

In summary, bears, wolves, and red foxes show similar phylogeographic structuring and evidence of large-scale gene flow, but also of recently reduced levels of connectivity and local adaptations in some regions.

Outlook

Fruitful future research will be a refined sampling in several geographic regions, especially in northern Africa, Asia, and Eastern Europe. Because some southern populations appear not to have contributed to large-scale postglacial range expansions, those seem particularly likely to harbor previously undetected genetic variation. High sequence variability in some mitochondrial genomic regions can provide enough information to detect phylogeographic events, especially the hypervariable 5′ end of the mammalian mtDNA control region. However, analysis of a larger mtDNA fragment or of the whole mitochondrial genome can reveal additional phylogeographic structure, in particular among recently diverged lineages [58-63]. Such a larger fragment could help overcome some of the topological uncertainties present in our dataset.

Microsatellite markers have been used to study the fine-scale population structure in geographically restricted red fox populations [24,52-54]. Frati et al. [48] used allozymes and cytochrome b sequences to compare genetic variability among some European red fox populations. To date, most studies on the large-scale population structuring and phylogeography of red foxes have utilized mtDNA sequences [21-27,29]. As a maternally inherited molecule with a high mutation rate, compared to the nuclear genome, and fast coalescence (due to lower effective population size than autosomal loci), mtDNA has been used to resolve phylogeographic structures in many taxa. However, as male-mediated gene flow cannot be detected from mtDNA, more complete inferences of the phylogeographic history of a species should include biparentally or paternally inherited markers (e.g., [64,65]), especially in species like the red fox where males disperse more than females [20].

Therefore, the application of SNP chips, modern high-throughput sequencing techniques [63], and/or the establishment of new (nuclear) markers in a synthesis with phenotypic features and ecological adaptations (e.g., [12,66]) will lead to a deeper understanding of the phylogeographic history and adaptations of this widespread generalist.

Methods

Samples and DNA extraction

In total, 52 red fox samples and one arctic fox sample were used in this study (Figure  1, Additional file 1). 33 red fox muscle and skin samples were collected in Germany, Finland, Poland, and central Siberia. The arctic fox muscle sample was collected in central Siberia (Additional file  1). All samples originated from dead wildlife legally hunted during hunting season according to local law, and the animals were not killed specifically for this study. No ethical approval or permit for animal experimentation was required. Additionally, 16 hair and 3 fecal samples were collected in Germany.

Total DNA was extracted from muscle and skin samples using a standard salt extraction protocol (modified from the Puregene™ DNA extraction kit). We extracted DNA from hair using the QIAamp DNA Investigator kit (Qiagen, Hilden) protocols for hair and the QIAamp DNA Stool Mini Kit for fecal samples following the manufacturers’ instructions.

Amplification and sequencing

A 449 bp fragment (excluding primers) from the 5′ end of the mitochondrial control region was amplified using the primers Vv.CRS1F 5′-CCCCAAGACTCAAGGAAGAGGCAC and Vv.CRS1R 5′-ACACCACAGTTATGTGTGATCATGGGC. These primers were newly designed based on an alignment of published Vulpes mitochondrial sequence [GenBank:GQ374180, AF098155, EU935091 (unpublished); GenBank:AM181037, NC_008434 [67]; GenBank:AB292765, AB292754, AB292741 [21]; GenBank:D83639 [68]; GenBank:AJ585358 [69]]. The forward primer is located at the end of the mitochondrial tRNA-Thr and the beginning of the tRNA-Pro genes, and the reverse primer is located 397 nucleotides into in the control region. The amplified region encompasses the hypervariable region targeted in previous red fox control region studies, allowing direct comparison of the data. PCR reactions were carried out in 15 μl volumes containing approximately 15 ng of genomic DNA, 0.27 μM of each primer, 0.16 μg/μl BSA (New England Biolabs, Ipswich, MA, USA), and 0.8 × of VWR Taq DNA Polymerase Master Mix containing a final concentration of 1.6 mM MgCl2 (VWR, Darmstadt, Germany). PCR was performed on a Unocycler (VWR, Darmstadt, Germany) using the following thermal profile: 3 min at 95°C prior to 40 cycles of 30 s at 94°C, 25 s at 59°C, and 1 min 15 s at 72°C; followed by an extension step of 10 min at 72°C. PCR products were detected using standard 1.5% agarose gel electrophoresis, and cycle sequenced with BigDye 3.1 chemistry (Applied Biosystems, Foster City, CA, USA) according to the manufacturer’s recommendation using 1/12th of the reaction mix, with 0.16 μl of BigDye in a 10 μl total volume reaction. Detected PCR products were run on an ABI 3100 instrument (Applied Biosystems). Electropherograms were assembled, checked manually, and sequences were aligned using Geneious 5.4 (Biomatters).

Sequence analyses

A total of 677 previously published and publicly available sequences from wild red fox populations were collated from GenBank for the control region [21-25,27,28] (Figure  1, Additional file 1). Aubry et al. [23], Sacks et al. [24], and Statham et al. [27] obtained most sequences from museum specimens (1850–1991), the 35 sequences from Teacher et al. [25] were from the Late Pleistocene and early Holocene. All other sequences were of recent origin (>1989) [21,22,24,27,28]. The 52 newly obtained red fox sequences (from samples collected since 2009) were added to a final alignment containing in total 729 red foxes, with a length of 335 bp. Ancient DNA samples from Teacher et al. [25] and one modern sequence from Valière et al. [28] were shorter (179–268 bp), and therefore omitted in network analyses (see below). A comprehensive list of all analyzed sequences, including geographic origins (see also Figure  1), GenBank accession numbers, and referenced study is provided in Additional file 1.

Data analyses

We determined parameters of within-population variability for all samples and major lineages using Arlequin 3.5 [70], calculating the number of segregating sites and haplotypes, haplotype diversity (Hd), nucleotide diversity (π), and Fu’s FS[71], an indicator of population expansion when it is negative and significant (p ≤ 0.02) [35,70,71] (Table  1).

Mismatch distributions under the sudden expansion model [72] were modeled and investigated in Arlequin 3.5 [70]. The sudden expansion model assumes population growth from a population at equilibrium with θ = θ0 to a new size with θ = θ1 within τ units of mutational time, with τ = 2*u*t (u = substitution rate per lineage for the entire DNA fragment, and t = number of generations since the expansion; see [73]). Time since expansion was thus calculated by dividing the estimate of τ by the product of: 335 (sequence length in base pairs) and the divergence rate (twice the per-lineage substitution rate; see [74]) in percent per year. Generation time for red foxes was assumed to be one year [16,75]. To test goodness-of-fit of the observed mismatch distribution to that expected under the sudden population expansion model, the sum of squared deviations [76] was computed in Arlequin with 10,000 replicates. To calculate time since expansion for those lineages that did not deviate from sudden expansion (p > 0.05), we assumed a per-lineage substitution rate of: (i) 7.1% substitutions per lineage per million years (or u = 2.379*10-5), a rate previously used for the arctic fox [2], and (ii) 33.2% substitutions per lineage per million years (u = 1.112*10-4), the mean rate estimated from two independent simulations in BEAST (Scenario 1; details see below) (Table  3).

A median-joining network of nucleotide sequences (n = 693) was constructed using the software Network 4.6.1.0 [77] (Figure  3). Due to their shorter sequence length and thus missing data, all 35 ancient samples from Teacher et al. [25] and one modern sample from Valière et al. [28] were excluded.

To identify the model of sequence evolution that best fit the data, we used jmodeltest 0.0.1 [78], which suggested the TN93 + G + I model of sequence evolution. A phylogeny and divergence time estimates for different lineages were obtained from simulations in BEAST v1.7.5 [79]. For computational reasons, a maximum number of eight individuals was included per haplotype. BEAST was set to run for 1 billion generations, sampling every 10,000th generation. Convergence was checked in Tracer v1.5. Two runs with identical settings were combined before resampling ca. 20,000 trees, both using LogCombiner v1.7.5 (without setting a burn-in). A maximum clade credibility tree was constructed using TreeAnnotator with a burn-in of 10%. Besides the constant population size model, we also evaluated the exponential population growth model implemented in BEAST.

A combined approach was used for dating of phylogeographic events, utilizing the arctic fox as exterior calibration point, plus tip dates based on published ancient DNA red fox sequences [25]. This approach accounts for major discrepancies from the possible time dependency of the molecular clock on recent evolutionary time scales [32]. The fossil tip calibrations are closer to the phylogeographic time frames of interest than the exterior calibration point, avoiding possible rate shifts. To further validate the robustness of our dating methodology, we performed a simulation using only tip dates, without setting a prior for the root height of the tree.

The divergence time to the arctic fox was used as exterior calibration point, because to our knowledge this is the closest relative to the red fox with an available divergence time estimate.

We tested three scenarios in BEAST, all three using the same tip dates, but varying root heights as exterior calibration points (Table  2). Scenario 1: the uniform prior for the root height was set to 1.75-4 million years (My), according to the 95% credibility interval in Perini et al. [80] for the divergence time between red and arctic fox (Figure  2). Scenario 2: we tested a uniform prior that spanned a period of 5.1-6 My, according to the 95% credibility interval for another divergence time estimate by Nyakatura and Bininda-Emonds [81]. Scenario 3: a very recent divergence time scenario was tested, based on the first appearance of the red fox in the fossil record (0.5-1 Mya) [30,31], setting the minimum age of the root height to 0.5 My and using a lognormal prior. As the speciation event very likely happened earlier than 0.5 Mya [80,81], the 95% interval of the lognormal prior included a period of up to 5.9 My.

A phylogenetic analysis of haplotype data was conducted in MrBayes 3.2 [82] (Additional file 3). We used jmodeltest 0.0.1 [78] to find the model of evolution that best fits the data (HKY + G + I). The analysis was run for 15 million MCMC generations sampling every 2,000th generation, with a burnin of 25%. We used four heated chains and confirmed convergence using the potential scale reduction factor [83], a convergence diagnostic implemented in MrBayes, which approached 1.000 for all parameters.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

VEK performed laboratory analyses, analyzed the data, and wrote the manuscript. AJ helped conceive the study and contributed reagents. NL, NS, AAS, TH, and CN contributed samples. KS contributed samples and extracted non-invasive samples. NL, AJ, NS, AAS, TH, KS, and CN participated in revising the manuscript. FH conceived the study, analyzed the data, and wrote the manuscript. All authors read and approved the final manuscript.

Acknowledgements

We thank Vasiliy Sokolov, Natalya Sokolova, Kamila Plis, Filip Zięba, Jochen Becker, and Tomas Häggvik for providing us with samples, and Anna Kaempffe and Jessica Appel for assistance in the lab. Sampling was partly funded by the Research Council of Norway within the IPY project "Arctic predators" and by the Program of UD RAS, project No 12-4-7-022-Arctic. We acknowledge the contribution of the Department of Environment/Government of Nunavut. We thank Frank Zachos and three anonymous reviewers for helpful comments on a previous version of the manuscript. The study was supported by Hesse’s funding program LOEWE (Landes-Offensive zur Entwicklung Wissenschaftlich-ökonomischer Exzellenz).

References

  1. Paetkau D, Amstrup SC, Born EW, Calvert W, Derocher AE, Garner GW, Messier F, Stirling I, Taylor MK, Wiig Ø, Strobeck C: Genetic structure of the world’s polar bear populations.

    Mol Ecol 1999, 8:1571-1584. PubMed Abstract | Publisher Full Text OpenURL

  2. Dalén L, Fuglei E, Hersteinsson P, Kapel CMO, Roth JD, Samelius G, Tannerfeldt M, Angerbjörn A: Population history and genetic structure of a circumpolar species: the arctic fox.

    Biol J Linn Soc 2005, 84:79-89. OpenURL

  3. Wooding S, Ward R: Phylogeography and pleistocene evolution in the North American black bear.

    Mol Biol Evol 1997, 14:1096-1105. PubMed Abstract | Publisher Full Text OpenURL

  4. Seddon JM, Santucci F, Reeve NJ, Hewitt GM: DNA footprints of European hedgehogs, Erinaceus europaeus and E. concolor: Pleistocene refugia, postglacial expansion and colonization routes.

    Mol Ecol 2001, 10:2187-2198. PubMed Abstract | Publisher Full Text OpenURL

  5. Niedziałkowska M, Jędrzejewska B, Honnen A-C, Otto T, Sidorovich VE, Perzanowski K, Skog A, Hartl GB, Borowik T, Bunevich AN, Lang J, Zachos FE: Molecular biogeography of red deer Cervus elaphus from eastern Europe: insights from mitochondrial DNA sequences.

    Acta Theriol 2011, 56:1-12. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  6. Taberlet P, Fumagalli L, Wust-Saucy A-G, Cosson J-F: Comparative phylogeography and postglacial colonization routes in Europe.

    Mol Ecol 1998, 7:453-464. PubMed Abstract | Publisher Full Text OpenURL

  7. Hewitt G: The genetic legacy of the Quaternary ice ages.

    Nature 2000, 405:907-913. PubMed Abstract | Publisher Full Text OpenURL

  8. Vilà C, Amorim IR, Leonard JA, Posada D, Castroviejo J, Petrucci-Fonseca F, Crandall KA, Ellegren H, Wayne RK: Mitochondrial DNA phylogeography and population history of the grey wolf Canis lupus.

    Mol Ecol 1999, 8:2089-2103. PubMed Abstract | Publisher Full Text OpenURL

  9. Sharma DK, Maldonado JE, Jhala YV, Fleischer RC: Ancient wolf lineages in India.

    P R Soc B 2004, 271:S1-S4. Publisher Full Text OpenURL

  10. Valdiosera CE, García N, Anderung C, Dalén L, Crégut-Bonnoure E, Kahlke R-D, Stiller M, Brandström M, Thomas MG, Arsuaga JL, Götherström A, Barnes I: Staying out in the cold: glacial refugia and mitochondrial DNA phylogeography in ancient European brown bears.

    Mol Ecol 2007, 16:5140-5148. PubMed Abstract | Publisher Full Text OpenURL

  11. Calvignac S, Hughes S, Hänni C: Genetic diversity of endangered brown bear (Ursus arctos) populations at the crossroads of Europe, Asia and Africa.

    Divers Distrib 2009, 15:742-750. Publisher Full Text OpenURL

  12. Muñoz-Fuentes V, Darimont CT, Wayne RK, Paquet PC, Leonard JA: Ecological factors drive differentiation in wolves from British Columbia.

    J Biogeogr 2009, 36:1516-1531. Publisher Full Text OpenURL

  13. Davison J, Ho SYW, Bray SC, Korsten M, Tammeleht E, Hindrikson M, Østbye K, Østbye E, Lauritzen S-E, Austin J, Cooper A, Saarma U: Late-Quaternary biogeographic scenarios for the brown bear (Ursus arctos), a wild mammal model species.

    Q Sci Rev 2011, 30:418-430. Publisher Full Text OpenURL

  14. Edwards CJ, Suchard MA, Lemey P, Welch JJ, Barnes I, Fulton TL, Barnett R, O’Connell TC, Coxon P, Monaghan N, Valdiosera CE, Lorenzen ED, Willerslev E, Baryshnikov GF, Rambaut A, Thomas MG, Bradley DG, Shapiro B: Ancient Hybridization and an Irish Origin for the Modern Polar Bear Matriline.

    Curr Biol 2011, 21:1251-1258. PubMed Abstract | Publisher Full Text OpenURL

  15. Weckworth BV, Dawson NG, Talbot SL, Flamme MJ, Cook JA: Going Coastal: Shared Evolutionary History between Coastal British Columbia and Southeast Alaska Wolves (Canis lupus).

    PLoS One 2011, 6:e19582. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  16. Larivière S, Pasitschniak-Arts M: Vulpes vulpes.

    Mamm Species 1996, 537:1-11. OpenURL

  17. Sommer R, Benecke N: Late-Pleistocene and early Holocene history of the canid fauna of Europe (Canidae).

    Mamm Biol 2005, 70:227-241. OpenURL

  18. Voigt DR: Red fox. In Wild furbearer management and conservation in North America. Edited by Novák M, Baker JA, Obbard ME, Malloch B. Toronto: Ontario Ministry of Natural Resources; 1987:379-382. OpenURL

  19. Goszczynski J: Population dynamics of the red fox in central Poland.

    Acta Theriol 1989, 34:141-154. OpenURL

  20. Allen SH, Sargeant AB: Dispersal Patterns of Red Foxes Relative to Population Density.

    J Wildl Manage 1993, 57:526-533. Publisher Full Text OpenURL

  21. Inoue T, Nonaka N, Mizuno A, Morishima Y, Sato H, Katakura K, Oku Y: Mitochondrial DNA Phylogeography of the Red Fox (Vulpes vulpes) in Northern Japan.

    Zool Sci 2007, 24:1178-1186. PubMed Abstract | Publisher Full Text OpenURL

  22. Kirschning J, Zachos FE, Cirovic D, Radovic IT, Hmwe SS, Hartl GB: Population Genetic Analysis of Serbian Red Foxes (Vulpes vulpes) by Means of Mitochondrial Control Region Sequences.

    Biochem Genet 2007, 45:409-420. PubMed Abstract | Publisher Full Text OpenURL

  23. Aubry KB, Statham MJ, Sacks BN, Perrine JD, Wisely SM: Phylogeography of the North American red fox: vicariance in Pleistocene forest refugia.

    Mol Ecol 2009, 18:2668-2686. PubMed Abstract | Publisher Full Text OpenURL

  24. Sacks BN, Statham MJ, Perrine JD, Wisely SM, Aubry KB: North American montane red foxes: expansion, fragmentation, and the origin of the Sacramento Valley red fox.

    Conserv Genet 2010, 11:1523-1539. Publisher Full Text OpenURL

  25. Teacher AG, Thomas JA, Barnes I: Modern and ancient red fox (Vulpes vulpes) in Europe show an unusual lack of geographical and temporal structuring, and differing responses within the carnivores to historical climatic change.

    BMC Evol Biol 2011, 11:214. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  26. Edwards CJ, Soulsbury CD, Statham MJ, Ho SYW, Wall D, Dolf G, Iossa G, Baker PJ, Harris S, Sacks BN, Bradley DG: Temporal genetic variation of the red fox, Vulpes vulpes, across western Europe and the British Isles.

    Q Sci Rev 2012, 57:95-104. OpenURL

  27. Statham MJ, Sacks BN, Aubry KB, Perrine JD, Wisely SM: The origin of recently established red fox populations in the United States: translocations or natural range expansions?

    J Mammal 2012, 93:52-65. Publisher Full Text OpenURL

  28. Valière N, Fumagalli L, Gielly L, Miquel C, Lequette B, Poulle M-L, Weber J-M, Arlettaz R, Taberlet P: Long-distance wolf recolonization of France and Switzerland inferred from non-invasive genetic sampling over a period of 10 years.

    Anim Conserv 2003, 6:83-92. Publisher Full Text OpenURL

  29. Yu J-N, Han S-H, Kim B-H, Kryukov AP, Kim S, Lee B-Y, Kwak M: Insights into Korean Red Fox (Vulpes vulpes) Based on Mitochondrial Cytochrome b Sequence Variation in East Asia.

    Zool Sci 2012, 29:753-760. PubMed Abstract | Publisher Full Text OpenURL

  30. Kurtén B: Pleistocene Mammals of Europe. Chicago, Illinois: Aldine Transaction; 1968. OpenURL

  31. Wayne RK, Van Valkenburgh B, O’Brien SJ: Molecular distance and divergence time in carnivores and primates.

    Mol Biol Evol 1991, 8:297-319. PubMed Abstract | Publisher Full Text OpenURL

  32. Ho SYW, Phillips MJ, Cooper A, Drummond AJ: Time Dependency of Molecular Rate Estimates and Systematic Overestimation of Recent Divergence Times.

    Mol Biol Evol 2005, 22:1561-1568. PubMed Abstract | Publisher Full Text OpenURL

  33. Saarma U, Ho SYW, Pybus OG, Kaljuste M, Tumanov IL, Kojola I, Vorobiev AA, Markov NI, Saveljev AP, Valdmann H, Lyapunova EA, Abramov AV, Männil P, Korsten M, Vulla E, Pazetnov SV, Pazetnov VS, Putchkovskiy SV, Rõkov AM: Mitogenetic structure of brown bears (Ursus arctos L.) in northeastern Europe and a new time frame for the formation of European brown bear lineages.

    Mol Ecol 2007, 16:401-413. PubMed Abstract | Publisher Full Text OpenURL

  34. Leonard JA, Vilà C, Wayne RK: Legacy lost: genetic variability and population size of extirpated US grey wolves (Canis lupus).

    Mol Ecol 2005, 14:9-17. PubMed Abstract | Publisher Full Text OpenURL

  35. Ramos-Onsins SE, Rozas J: Statistical Properties of New Neutrality Tests Against Population Growth.

    Mol Biol Evol 2002, 19:2092-2100. PubMed Abstract | Publisher Full Text OpenURL

  36. Pilot M, Jedrzejewski W, Branicki W, Sidorovich VE, Jedrzejewska B, Stachura K, Funk SM: Ecological factors influence population genetic structure of European grey wolves.

    Mol Ecol 2006, 15:4533-4553. PubMed Abstract | Publisher Full Text OpenURL

  37. Kurtén B, Anderson E: Pleistocene Mammals of North America. North America: Columbia University Press; 1980. OpenURL

  38. Hoffecker JF, Elias SA: The Human Ecology of Beringia. New York: Columbia University Press; 2007. OpenURL

  39. Péwé TL, Hopkins DM: Mammal Remains of Pre-Wisconsin Age in Alaska. California: University Press; 1967. OpenURL

  40. Péwé TL: Quaternary Geology of Alaska. Washington, DC: US Government Printing Office; 1975. [Geological Survey Professional Paper no 835] OpenURL

  41. Harington CR: Pleistocene vertebrates of the Yukon Territory. University of Alberta; 1977. [Dissertation] OpenURL

  42. FAUNMAP Working Group: FAUNMAP: a database documenting late Quaternary distributions of mammal species in the United States.

    Illinois State Museum Scientific Papers 1994, 25:1-690. OpenURL

  43. McKay BD: A new timeframe for the diversification of Japan’s mammals.

    J Biogeogr 2011, 39:1134-1143. OpenURL

  44. Blakiston T: Zoological indications of ancient connection of the Japan Islands with the continent.

    Trans Asiatic Soc Japan 1883, 11:126-140. OpenURL

  45. Ding Q-L, Wang C-C, Farina SE, Li H: Mapping Human Genetic Diversity on the Japanese Archipelago.

    Adv Anthropol 2011, 1:19-25. Publisher Full Text OpenURL

  46. Kawamura Y: Last glacial and Holocene land mammals of the Japanese islands: their fauna, extinction and immigration.

    Q Res 2007, 46:171-177. Publisher Full Text OpenURL

  47. Hersteinsson P, MacDonald DW: Interspecific Competition and the Geographical Distribution of Red and Arctic Foxes Vulpes vulpes and Alopex lagopus.

    Oikos 1992, 64:505-515. Publisher Full Text OpenURL

  48. Frati F, Hartl GB, Lovari S, Delibes M, Markov G: Quaternary radiation and genetic structure of the red fox Vulpes vulpes in the Mediterranean Basin, as revealed by allozymes and mitochondrial DNA.

    J Zool 1998, 245:43-51. Publisher Full Text OpenURL

  49. Tannerfeldt M, Elmhagen B, Angerbjörn A: Exclusion by interference competition? The relationship between red and arctic foxes.

    Oecologia 2002, 132:213-220. Publisher Full Text OpenURL

  50. Killengreen ST, Lecomte N, Ehrich D, Schott T, Yoccoz NG, Ims RA: The importance of marine vs. human-induced subsidies in the maintenance of an expanding mesocarnivore in the arctic tundra.

    J Anim Ecol 2011, 80:1049-1060. PubMed Abstract | Publisher Full Text OpenURL

  51. Rodnikova A, Ims RA, Sokolov AA, Skogstad G, Sokolov V, Shtro V, Fuglei E: Red fox takeover of arctic fox breeding den: an observation from Yamal Peninsula, Russia.

    Polar Biol 2011, 34:1609-1614. Publisher Full Text OpenURL

  52. Robinson NA, Marks CA: Genetic structure and dispersal of red foxes (Vulpes vulpes) in urban Melbourne.

    Aust J Zool 2001, 49:589-601. Publisher Full Text OpenURL

  53. Wandeler P, Funk SM, Largiadèr CR, Gloor S, Breitenmoser U: The city-fox phenomenon: genetic consequences of a recent colonization of urban habitat.

    Mol Ecol 2003, 12:647-656. PubMed Abstract | Publisher Full Text OpenURL

  54. Oishi T, Uraguchi K, Takahashi K, Masuda R: Population Structures of the Red Fox (Vulpes vulpes) on the Hokkaido Island, Japan, Revealed by Microsatellite Analysis.

    J Hered 2011, 102:38-46. PubMed Abstract | Publisher Full Text OpenURL

  55. Korsten M, Ho SYW, Davison J, Pähn B, Vulla E, Roht M, Tumanov IL, Kojola I, Andersone-Lilley Z, Ozolins J, Pilot M, Mertzanis Y, Giannakopoulos A, Vorobiev AA, Markov NI, Saveljev AP, Lyapunova EA, Abramov AV, Männil P, Valdmann H, Pazetnov SV, Pazetnov VS, Rõkov AM, Saarma U: Sudden expansion of a single brown bear maternal lineage across northern continental Eurasia after the last ice age: a general demographic model for mammals?

    Mol Ecol 2009, 18:1963-1979. PubMed Abstract | Publisher Full Text OpenURL

  56. Matsuhashi T, Masuda R, Mano T, Murata K, Aiurzaniin A: Phylogenetic Relationships among Worldwide Populations of the Brown Bear Ursus arctos.

    Zool Sci 2001, 18:1137-1143. Publisher Full Text OpenURL

  57. Pilot M, Jędrzejewski W, Sidorovich VE, Meier-Augenstein W, Hoelzel AR: Dietary Differentiation and the Evolution of Population Genetic Structure in a Highly Mobile Carnivore.

    PLoS One 2012, 7:e39341. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  58. Björnerfeldt S, Webster MT, Vilà C: Relaxation of selective constraint on dog mitochondrial DNA following domestication.

    Genome Res 2006, 16:990-994. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  59. Lindqvist C, Schuster SC, Sun Y, Talbot SL, Qi J, Ratan A, Tomsho LP, Kasson L, Zeyl E, Aars J, Miller W, Ingólfsson Ó, Bachmann L, Wiig Ø: Complete mitochondrial genome of a Pleistocene jawbone unveils the origin of polar bear.

    Proc Natl Acad Sci USA 2010, 107:5053-5057. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  60. Morin PA, Archer FI, Foote AD, Vilstrup J, Allen EE, Wade P, Durban J, Parsons K, Pitman R, Li L, Bouffard P, Nielsen SCA, Rasmussen M, Willerslev E, Gilbert MTP, Harkins T: Complete mitochondrial genome phylogeographic analysis of killer whales (Orcinus orca) indicates multiple species.

    Genome Res 2010, 20:908-916. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  61. Rutledge LY, Patterson BR, White BN: Analysis of Canis mitochondrial DNA demonstrates high concordance between the control region and ATPase genes.

    BMC Evol Biol 2010, 10:215. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  62. Keis M, Remm J, Ho SYW, Davison J, Tammeleht E, Tumanov IL, Saveljev AP, Männil P, Kojola I, Abramov AV, Margus T, Saarma U: Complete mitochondrial genomes and a novel spatial genetic method reveal cryptic phylogeographical structure and migration patterns among brown bears in north-western Eurasia.

    J Biogeogr 2013, 40:915-927. Publisher Full Text OpenURL

  63. McCormack JE, Hird SM, Zellmer AJ, Carstens BC, Brumfield RT: Applications of next-generation sequencing to phylogeography and phylogenetics.

    Mol Phylogenet Evol 2013, 66:526-538. PubMed Abstract | Publisher Full Text OpenURL

  64. Hellborg L, Gündüz İ, Jaarola M: Analysis of sex-linked sequences supports a new mammal species in Europe.

    Mol Ecol 2005, 14:2025-2031. PubMed Abstract | Publisher Full Text OpenURL

  65. Hailer F, Kutschera VE, Hallstrom BM, Klassert D, Fain SR, Leonard JA, Arnason U, Janke A: Nuclear Genomic Sequences Reveal that Polar Bears Are an Old and Distinct Bear Lineage.

    Science 2012, 336:344-347. PubMed Abstract | Publisher Full Text OpenURL

  66. Nadachowska-Brzyska K, Zieliński P, Radwan J, Babik W: Interspecific hybridization increases MHC class II diversity in two sister species of newts.

    Mol Ecol 2012, 21:887-906. PubMed Abstract | Publisher Full Text OpenURL

  67. Arnason U, Gullberg A, Janke A, Kullberg M, Lehman N, Petrov EA, Väinölä R: Pinniped phylogeny and a new hypothesis for their origin and dispersal.

    Mol Phylogenet Evol 2006, 41:345-354. PubMed Abstract | Publisher Full Text OpenURL

  68. Okumura N, Ishiguro N, Nakano M, Matsui A, Sahara M: Intra- and interbreed genetic variations of mitochondrial DNA major non-coding regions in Japanese native dog breeds (Canis familiaris).

    Anim Genet 1996, 27:397-405. PubMed Abstract OpenURL

  69. Statham MJ, Turner PD, O’Reilly C: Use of PCR amplification and restriction enzyme digestion of mitochondrial D-loop for identification of mustelids in Ireland.

    Ir Nat J 2005, 28:1-6. OpenURL

  70. Excoffier L, Lischer HEL: Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows.

    Mol Ecol Resour 2010, 10:564-567. PubMed Abstract | Publisher Full Text OpenURL

  71. Fu Y-X: Statistical Tests of Neutrality of Mutations Against Population Growth, Hitchhiking and Background Selection.

    Genetics 1997, 147:915-925. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  72. Rogers AR, Harpending H: Population growth makes waves in the distribution of pairwise genetic differences.

    Mol Biol Evol 1992, 9:552-569. PubMed Abstract | Publisher Full Text OpenURL

  73. Rogers AR: Genetic Evidence for a Pleistocene Population Explosion.

    Evolution 1995, 49:608-615. Publisher Full Text OpenURL

  74. Schenekar T, Weiss S: High rate of calculation errors in mismatch distribution analysis results in numerous false inferences of biological importance.

    Heredity 2011, 107:511-512. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  75. Lloyd HG: The Red Fox. London: Batsford; 1980. OpenURL

  76. Schneider S, Excoffier L: Estimation of Past Demographic Parameters From the Distribution of Pairwise Differences When the Mutation Rates Vary Among Sites: Application to Human Mitochondrial DNA.

    Genetics 1999, 152:1079-1089. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  77. Bandelt HJ, Forster P, Röhl A: Median-joining networks for inferring intraspecific phylogenies.

    Mol Biol Evol 1999, 16:37-48. PubMed Abstract | Publisher Full Text OpenURL

  78. Posada D: jModelTest: Phylogenetic Model Averaging.

    Mol Biol Evol 2008, 25:1253-1256. PubMed Abstract | Publisher Full Text OpenURL

  79. Drummond AJ, Suchard MA, Xie D, Rambaut A: Bayesian phylogenetics with BEAUti and the BEAST 1.7.

    Mol Biol Evol 2012, 29:1969-1973. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  80. Perini FA, Russo CAM, Schrago CG: The evolution of South American endemic canids: a history of rapid diversification and morphological parallelism.

    J Evol Biol 2010, 23:311-322. PubMed Abstract | Publisher Full Text OpenURL

  81. Nyakatura K, Bininda-Emonds OR: Updating the evolutionary history of Carnivora (Mammalia): a new species-level supertree complete with divergence time estimates.

    BMC Biol 2012, 10:12. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  82. Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, Larget B, Liu L, Suchard MA, Huelsenbeck JP: MrBayes 3.2: Efficient Bayesian Phylogenetic Inference and Model Choice Across a Large Model Space.

    Syst Biol 2012, 61:539-542. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  83. Gelman A, Rubin DB: Inference from Iterative Simulation Using Multiple Sequences.

    Stat Sci 1992, 7:457-472. Publisher Full Text OpenURL

  84. IUCN (International Union for Conservation of Nature): Vulpes vulpes.

    IUCN 2011 IUCN Red List of Threatened Species Version 20121 2008. OpenURL