Across the Arctic: Mitogenomic Phylogeny of Arctic Foxes (Vulpes lagopus) Reveals Several New Matrilines and Illuminates the Colonization History of the Icelandic Population

TL;DR

This study uses mitochondrial DNA to trace the evolutionary history and colonization of Arctic foxes, including those in Iceland.

Contribution

The study identifies seven new haplogroups and provides insights into the colonization history of Icelandic Arctic foxes.

Findings

Seven distinct haplogroups were identified, diverging as far back as 171 thousand years ago.

Three haplogroups were found in Iceland, requiring at least four unrelated founding females.

Icelandic matrilineal diversity suggests multiple colonization sources but lacks a unique geographic origin.

Abstract

Background/Objectives: Arctic foxes (Vulpes lagopus) exemplify the vulnerability of Arctic species to global warming and anthropogenic impacts, including habitat loss, interspecific competition with temperate species, pollution (chemical and biological), and declining prey abundance. Despite their ecological importance, the evolutionary and demographic history of the species is still incompletely understood, and the colonization history of isolated island populations, such as the one on Iceland, remains unresolved. Methods: We analyzed 80 mitochondrial genomes from across the Holarctic, including 22 Icelandic individuals. We combined phylogenetic reconstruction, coalescence-dating, haplotype network analysis, and diversity metrics to infer matrilineal relationships and colonization history. Results: Seven distinct haplogroups (Hg.1–Hg.7) were identified, which diverged ≥65 thousand years ago (kya). Two haplogroups were broadly distributed across Fennoscandia, Russia, Iceland, and Canada, while others were region-specific: two in eastern Russia (respectively diverging ~171 kya and ~89 kya), one in central Russia (~66 kya), and two in Iceland (~95 kya and ~66 kya). Three haplogroups were detected in Iceland, and at least four unrelated founding females are required to explain the current matrilineal diversity. One haplogroup contained sufficient representatives for molecular dating, yielding a minimum colonization age of ~5600 years, assuming in situ diversification. Observed matrilineal diversity in Iceland does not uniquely identify a single geographic source. Conclusions: Arctic foxes’ distribution and diversity reflect repeated cycles of isolation and expansion as circumpolar environments shifted. Broader sampling across the Nearctic is critical to clarify the timing, sources, and routes of Iceland’s colonization, as Nearctic sampling was limited to a single Canadian mitogenome.

Figures

Full text

1. Introduction

The Arctic is one of the most rapidly changing ecosystems in the world [1]. It is defined by cold conditions, permafrost, and strong seasonality in the availability of light and biological resources. The region is an early indicator of global warming and its effects on biodiversity and ecosystem structure [2]. Because of its harshness, the Arctic supports a highly specialized flora and fauna. However, due to their adaptation to the extreme conditions of the Arctic, these species are particularly susceptible to anthropogenic pressures (or anthropogenically driven niche changes) [1]. ‘Arctic amplification’ has led to significantly faster warming in the region compared to the global average [3], resulting in widespread disruption of habitat structure, species interactions, and migration corridors [1,4,5]. Understanding how Arctic species have historically responded to climate oscillations will be crucial for predicting their resilience to current and future warming [6,7,8]. With its circumpolar distribution and sensitivity to environmental change, the Arctic fox (Vulpes lagopus) has been chosen as a climate change flagship species for this biome [9].

The Arctic fox is a highly specialized carnivore adapted to the extreme environment of the Arctic tundra. Its distribution includes Alaska, Canada, Greenland, Iceland, Scandinavia, and Russia. It possesses a dense, insulating coat with seasonally changing colouration, a compact body, and the capacity to slow its metabolic rate to survive starvation periods [10]. Arctic foxes are opportunistic feeders, exhibiting behavioral flexibility in response to temporal and spatial changes in prey abundance [11]. They are highly mobile, covering up to 6000 km in a year [12], and depend on sea ice for long-distance dispersal [7,13,14].

Climate change and anthropogenic activity pose major threats to this species. Interspecific competition with red foxes for denning sites and resources, climate-driven habitat loss, chemical and biological pollution, reduced dispersal ability due to declining sea ice, and declining prey abundance threaten Arctic foxes in parts of their range [13,14,15,16,17,18,19].

In Iceland, the Arctic fox population is naturally isolated from other populations by the oceanic barrier. The island’s landscape can be divided into two broad habitat types: coastal and inland [20]. ‘Coastal foxes’ predominantly inhabit the north-western and eastern regions [21], feeding on sea birds, fish- and seal carcasses, marine invertebrates, seaweeds and berries. These regions of the island offer aggregated resources and high productivity, supporting a high density of dens and strong territoriality among individuals [22]. In contrast, the diet of ‘inland foxes’ consists of ptarmigans (Lagopus muta), snow buntings (Plectrophenax nivalis), various migratory geese and waders, and wood mice (Apodemus sylvaticus) [20,23]. Due to relatively evenly distributed resources, inland foxes have larger territories and move longer distances to find food, while their populations are not strongly tied to cyclic prey dynamics like their continental counterparts, which are dependent on lemming population cycles [20,21]. Although Arctic foxes have been legally protected in Iceland since 1994 (Protection and Hunting of Wild Species Act [no. 64/1994]), historical foxhunting for livestock protection is still allowed. Following the population bottleneck in the late 1970s, the Arctic fox population size recovered at the same time as migrating bird populations expanded, in concurrence with rising temperature; despite this, Icelandic Arctic foxes exhibit low genetic diversity relative to other populations [21,24,25].

The history of how and when Arctic foxes colonized Iceland remains uncertain. A mitochondrial DNA study of 191 Arctic foxes from throughout the species range [25] indicated a global population expansion after the last interglacial maximum, ~118 thousand years ago (kya), followed by broad distribution during the subsequent glacial period when geographical barriers were minimal. Formerly glaciated regions, including Iceland, Greenland, Svalbard and Fennoscandia, were colonized after the Last Glacial Maximum (22 kya), and little phylogeographic structuring was detected across most of the species’ range, possibly due to the Arctic foxes’ strong dispersal ability [25]. Exceptions involved oceanic islands such as Iceland, where the lack of seasonal or permanent sea ice may have restricted gene flow and permitted divergence from other populations [13]. However, only a relatively short portion of the mitochondrial (mt) control region (292 nt) was available for these analyses [25]. Subsequent studies employing longer mtDNA sequences [26,27] did not include samples from Iceland, but revealed the existence of two major mtDNA clades in extant samples, albeit without geographic structuring, and evidence for additional lineages from archival samples.

Analysing complete mitochondrial genomes (‘mitogenomes’) instead of short mtDNA fragments considerably improves data resolution, enabling the reconstruction of colonization dynamics and the resolution of fine-scale phylogeographic patterns. As demonstrated in studies on the biogeography of other mammals, such as the brown bear [28] and the wolf [29,30], complete mitogenomes enable the detection of subtle population dynamics and allow the reconstruction of species’ demographic history that remains obscured in short DNA fragment-based analyses. Furthermore, higher data resolution is especially important for species in which low mtDNA diversity is expected [28]. Reconstructing the Arctic fox’s demographic and phylogeographic history using complete mitogenomes thus provides an opportunity to improve our understanding of how Arctic mammals persist and diversify under shifting environmental conditions.

This study aims to understand the colonization history of Iceland by Arctic foxes by analyzing complete mitogenomes from across a wide range of the Arctic foxes’ distribution. We identify matrilines, reconstruct their phylogenetic relationships and divergence times, and infer the number and timing of Iceland’s colonization events. This framework allows us to determine the source and origin of Iceland’s colonization by Arctic foxes and whether mitochondrial diversity in Iceland is structured by geography and/or ecotype.

2. Materials and Methods

2.1. Sampling and Data Sources

We collected organ samples (liver, kidney, muscle) from hunted Arctic foxes in the 2011 and 2012 hunting seasons from various inland and coastal regions of Iceland [31,32]. Legal hunting is managed by the Environmental Agency under the Wildlife Act (64/1994) and the Animal Welfare Act (55/2013). The scientific use of the hunting bag was organized by the Natural Science Institute of Iceland, courtesy of the Icelandic Ministry for the Environment and Natural Resources. Carcasses were stored at −20 °C until necropsy was performed. We used liver samples from 22 individuals for the genome analysis, with approximately a 1:1 ratio between the two ecotypes (Table 1).

In addition to the Icelandic samples, we generated complete mitochondrial genomes from WGS data for 58 samples: 40 from Sweden, four from Norway, 13 from Russia, and one from Canada [26,33,34]. Thus, in total, we analyzed mitogenomes from 80 Arctic fox samples (Figure 1; Supplementary Table S1). Additionally, 12 mitogenomes (Supplementary Table S2) from other canid species, mostly from the genus Vulpes, were downloaded from the NCBI database and included in the interspecific phylogenomic analysis.

2.2. DNA Extraction, Library Preparation, Sequencing, and Read Processing

DNA of Icelandic arctic foxes was extracted using a commercial DNA extraction kit (GENIAL GmbH, Troisdorf, Germany) following the manufacturer’s instructions. PCR-free libraries were sequenced on the Illumina NovaSeq 6000 platform (BMK, Beijing, China) with 150-bp paired-end reads.

Raw sequencing reads of all 80 samples were adapter-clipped, quality-trimmed and filtered using fastp v0.24.0 [36], using --qualified_quality_phred 15 --unqualified_percent_limit 40 --average_qual 30 -l 50. Pre- and post-trimming read quality was assessed using FastQC v0.12.1 [37] to verify overall high sequence quality and the effectiveness of adapter and quality trimming. Clipped, trimmed, and filtered reads were then mapped against a V. lagopus reference mitogenome (GenBank acc. no. NC_026529.1) using bwa-mem2 v2.2.1 [38]. SAM files were converted to BAM files using sambamba v1.0.1 [39], with the option -F “not unmapped and not mate_is_unmapped and mapping_quality ≥ 30 and sequence_length ≥ 50”. Sorting and duplicate removal were carried out with sambamba sort and sambamba markdup. Additionally, soft-clipped reads were removed using samtools view (v1.22.1) and awk.

2.3. Consensus Generation and Alignment

Consensus sequences for all 80 samples were generated using Geneious R8.1.9 [40]. Sites were N-masked if their coverage was below 20× or when they did not conform to a 75% majority rule. The 80 consensus sequences obtained were aligned using MAFFT v7.505 [41], and ambiguous or poorly aligned positions were manually curated. A 242-nt-long repetitive portion in the control region (positions 16095 to 16336 in the reference) was removed because it could not be reliably resolved with short-read data, resulting in a final alignment of 16414, not without ambiguities (i.e., no Ns).

2.4. Summary Statistics

Summary statistics for the entire dataset and subsets thereof were estimated using DnaSP v6.12.03 [42], including the number of segregating sites (S), number of haplotypes (H), haplotype diversity (Hd), nucleotide diversity (Pi), and average number of nucleotide differences (k). Resulting values were validated using 10,000 permutations. Diversity indices were not estimated for populations represented by single samples (Canada, Kola peninsula).

2.5. Phylogenetic Analyses

To calibrate the evolutionary timescale for the Arctic fox, we reconstructed a phylogeny using complete mitochondrial genomes from members of the tribe Vulpini and three other canid taxa (wolf Canis lupus, raccoon dog Nyctereutes procyonoides, and bat-eared fox Otocyon megalotis), along with the two most divergent Arctic fox sequences (ERR10447078; ERR10447093). This alignment included kit fox (V. macrotis, PV090784.1), corsac fox (V. corsac, NC_023958.1), Rüppell’s fox (V. rueppellii, NC_087722.1), Tibetan fox (V. ferrilata, NC_027935.1), swift fox (V. velox, PV090785.1), red fox (V. vulpes, NC_008434.1), Cape fox (V. chama, NC_070063.1), Blanford’s fox (V. cana, PV090783.1), fennec fox (V. zerda, NC_023459.1), wolf (C. lupus, NC_008092.1), raccoon dog (N. procyonoides, NC_013700.1), and bat-eared fox (O. megalotis, NC_036369.1). As some of these sequences lacked the D-loop, this region was removed from this analysis, giving a total alignment length of 15,461 nt.

Partitions and their substitution models were selected with PartitionFinder2 [43] using the greedy algorithm and the Bayesian Information Criterion (BIC), which favors parsimonious yet well-fitting models. PartitionFinder2 identified eight subsets and their substitution models (see Supplementary Table S3). We used BEAST v2.7.7 [44] to reconstruct an interspecific, time-calibrated phylogeny under a relaxed lognormal clock and a birth-death tree prior, with an initial value of 0.01 substitutions per site per million years. Divergence times were constrained using uniform MRCA priors based on the canid fossil record, enforcing monophyly for all calibrated clades. The following calibrations (in million years ago, Mya) were applied: (1) the stem split between the tribes Canini and Vulpini (7–14 Mya; [45], (2) the MRCA of V. corsac and V. ferrilata (0.8–5 Mya; [46], (3) the MRCA of V. lagopus, V. macrotis, and V. velox (0.4–4 Mya, [47], and (4) the MRCA of V. rueppellii and V. vulpes (0.6–4 Mya; [48]. The two V. lagopus sequences were constrained as monophyletic. All other priors and model parameters not explicitly mentioned were left at their default values. Convergence and mixing of MCMC parameters were verified in Tracer v1.7.2 [49]. Eight independent runs of 40 million MCMC generations were combined using LogCombiner v2.7 and summarized with TreeAnnotator v2.7 (both included in the BEAST v2.7.7 package [44], discarding the first 75% as burn-in. A maximum clade credibility (MCC) tree was generated in TreeAnnotator v2.7.7 and visualized in FigTree v1.4.4 [50]; Supplementary Figure S1).

Building on the divergence estimates obtained from the interspecific analysis, a time-calibrated intraspecific phylogeny of Arctic foxes was reconstructed using the complete mitogenome dataset, including the control region (total length: 16,414 nt). This dataset includes predominantly Palearctic samples, with Nearctic representation limited to a single Canadian individual. As this restricts our ability to detect Nearctic-specific lineages, we note that coalescent-based inferences may be biased towards the demographic history of Palearctic populations. PartitionFinder2 identified five subsets and their substitution models (see Supplementary Table S4). Analyses were conducted in BEAST v2.7.7 under a strict clock model and a Coalescent Bayesian Skyline tree prior. A strict molecular clock was used due to the expected low variation in substitution rate at shallow evolutionary timescales within a single species. The clock rate prior was parameterized as lognormal (mean = −4.1605, stdev = 0.156), determined in the interspecific analysis (mean = 0.0159, stdev = 2.57 × 10^−3^). We applied a lognormal MRCA prior to the root, with parameters: mean = −1.4422, and stdev = 0.3579 (median 0.2498 Mya). Eight independent MCMC chains of 30 million generations were run, and convergence was confirmed by effective sample sizes (ESS > 200) in Tracer v1.7.2. The first 75% of each run were discarded as burn-in, and the post burn-in trees were combined using LogCombiner v2.7. A maximum clade credibility (MCC) tree was summarized with TreeAnnotator v2.7.7 using median node heights and visualized in FigTree v1.4.4. All other priors and operator settings not explicitly mentioned were retained at their default values.

Maximum-likelihood (ML) phylogenetic inferences were conducted in RAxML v8.2.12 [51] for both the inter- and intra-specific datasets. PartitionFinder2 was used to determine partitions and their substitution models for both datasets (Supplementary Tables S5 and S6) using the Bayesian Information Criterion (BIC) and a greedy search algorithm. For both datasets, rapid bootstrapping was performed jointly with the best-scoring ML tree search in two runs, using 100 bootstrap replicates each to evaluate node support (Supplementary Figures S2 and S3).

To further explore mitochondrial diversity and to identify distinct lineages within V. lagopus, a Median-Joining (MJ) network [52] was constructed in PopART v1.7 [53].

3. Results

We successfully recovered complete mitochondrial genomes for all 80 Arctic foxes analyzed (Supplementary Table S1), including the 22 samples from Iceland (Table 1). Phylogenetic analyses revealed the presence of seven distinct and well-supported mitochondrial lineages, irrespective of method employed (MJ-network, Figure 2; Bayesian and ML phylogenies, Figure 3, Supplementary Figures S3 and S4). Because previous mitochondrial DNA studies of Arctic foxes used different nomenclatures for the two previously described major lineages [26,27], and were based on substantially shorter mitochondrial sequences, we propose a revised naming scheme, designating the seven lineages as haplogroups 1–7 (abbreviated Hg.1–Hg.7). In the MJ network, haplogroups differ by at least 38 substitutions (Figure 2), and the dated Bayesian phylogeny suggests minimum divergence times of ≥65 thousand years ago (kya) from their closest sister lineages (Figure 3).

The root age for V. lagopus inferred from these 80 mitogenomes is 171 kya (95% HPD: 101–247), suggesting that the major mitochondrial lineages of V. lagopus originated during the Penultimate Glacial Period (MIS 6) in the Middle Pleistocene. This estimate is older than estimates of 118 kya in a species-wide study using shorter mitochondrial sequences (292 nt; [25], or of 87 kya in a study using longer sequences that included archival material (3954 nt, [26]. Our estimate is similar to that of 190 kya in another study employing archival material (3044 nt, [27]).

We note that our results reflect mostly Palearctic mitochondrial variation, with only a single mitogenome from the Nearctic (Canada). Below, we describe each haplogroup in phylogenetic and geographic context, including their relative ages, spatial distributions, and correspondence with major climatic phases.

Hg.1 comprises a single sequence from an eastern Russian sample (RUS_east_01) and represents the earliest diverging lineage in the dataset. It is separated from all other haplogroups by at least 93 substitutions and diverged at approximately 171 kya (95% HPD = 101–247 kya). This corresponds to the Penultimate Glacial Period (MIS 6), when Eurasian ice sheets expanded to their near-maximum extent, sea level dropped by ~120 m [56], and periglacial steppe stretched southward to ~48° N [57].

Hg.2 diverged from Hg.3–Hg.7 approximately 123 kya (95% HPD = 69–184 kya), and comprises two lineages that diverged ~47 kya (95% HPD = 23–73 kya). The first includes 28 Fennoscandian samples (27 Swedish and one western Russian). Twenty-four of the Swedish samples share a haplotype, while the remaining three have distinct but closely related sequences, differing by ≤two substitutions. The second lineage consists of one sequence from a central Russian sample (RUS_central_09), which is separated from the Fennoscandian cluster by ≥29 mutational steps. The divergence of Hg.2 from Hg.3–Hg.7 falls within the Last Interglacial (MIS 5e, 130–115 kya), just after the Eemian interglacial (~125 kya), a warm phase with global temperatures +2–4 °C above modern levels, when forest–tundra extended northward and continental ice volume was at a minimum [58]. The internal divergence ~47 kya (95% HPD = 23–73 kya) occurred during MIS 3 (60–27 kya), characterized by repeated cold stadials and interstadials, partial ice retreats, and the expansion of open tundra–steppe across Fennoscandia and northern Russia [57].

The phylogeny reveals a major bifurcation separating the clade containing Hg.3 and Hg.4 from that comprising Hg.5, Hg.6 and Hg.7 approximately 100 kya (95% HPD = 70–184 kya). The branching order here is not resolved with confidence (posterior probability < 0.8; ML bootstrap < 80), but the grouping of these haplogroups is consistent across Bayesian and maximum-likelihood analyses (Figure 3, Supplementary Figure S3). These divergences occurred after the interstadial in MIS 5c, when global mean temperature declined and steppe–tundra expanded again across northern Eurasia [57].

Hg.3 is the most diverse and geographically widespread haplogroup, encompassing samples from Fennoscandia (Sweden n = 13; Norway n = 4), central Russia (n = 5), eastern Russia (n = 1), Iceland (n = 5), and Canada (n = 1). Its divergence from Hg.4 occurs shortly after the divergence from Hg.5, Hg.6, and Hg.7, around 95 kya (95% HPD 53–141 kya), during the cooling climate of MIS 5b. Around 55 kya (95% HPD 30–81 kya) in MIS 3, Hg.3 exhibits pronounced internal diversification, but the branching order is not resolved with high confidence (posterior probability < 0.8; ML bootstrap < 80). This unresolved topology likely reflects short internode lengths (often associated with rapid diversification) and limited phylogenetic signal in our mitogenome dataset for resolving these closely spaced divergences. This period is marked by predominantly small oscillations of cold climate, expansion of continental ice sheets across northern Europe, exposure of continental shelves, and the spread of open, extensive, discontinuous steppe-tundra across northern Eurasia [59]. Iceland is represented by two haplotypes: one shared by four samples (ICE_inland_02, 03, 06, 07) and one (ICE_inland_05) that is closely related to an eastern Russian sample (RUS_east_03). Within Fennoscandia, three haplotypes were found in eight, five and four samples, respectively. Two of these differ by only three substitutions and cluster with two Russian sequences (RUS_central_03, 04). All seven Russian sequences in Hg.3 have distinct haplotypes, with RUS_central_01 and RUS_central_06 diverging from the rest of Hg.3~47 kya (95% HPD = 25–70 kya). The single Canadian sequence (CAN_01) is most closely related to RUS_central_02, both diverging from the rest of the haplogroup ~49 kya (95% HPD = 26–74 kya).

Hg.4 consists of only Icelandic samples (n = 5) that have three closely related haplotypes, differing by ≤two substitutions. The divergence from Hg.3 occurred ~95 kya (95% HPD = 53–141 kya), after the MIS 5c interstadial.

Hg.5 diverged from Hg.6 and Hg.7 ~89 kya (95% HPD = 49–133 kya) and is represented by a single eastern-Russian sample (RUS_east_02). In the Bayesian and maximum-likelihood phylogenies, this divergence occurred shortly after the split between [Hg.3 and Hg.4] and [Hg.5, Hg.6 and Hg.7], and its placement is not resolved with confidence (posterior probability < 0.8; ML bootstrap < 80). This period (MIS 5b) was characterized by progressive cooling, expansion of permafrost, and widespread development of arid steppe habitats, while northeastern Siberia remained largely unglaciated [60].

Hg.6 diverged from Hg.7 around 66 kya (95% HPD = 35–100 kya), within MIS 4 and before the sustained warming that marks the MIS 4/3 transition (~57 kya). This interval featured a major glacial expansion in northern Europe, lowered sea levels with exposed continental shelves, and largely open, discontinuous steppe–tundra across unglaciated northern Eurasia [57]. Hg.6 consists of only Icelandic samples (n = 12), representing eight closely related haplotypes separated by ≤two mutations. The coalescence time for this haplogroup is 5.6 kya (95% HPD = 1.4–11 kya); if this represents in situ diversification following colonization of Iceland, it provides a minimum age for occupancy by Arctic foxes.

Hg.7 consists of two closely related haplotypes (2 substitutions) from central Russian samples.

Across haplogroups, shallow internal splits (<10 mutational steps) occurred during the terminal Pleistocene and early Holocene (<15 kya; MIS 2–1 transition). These tip-level divergences coincide with the end of the LGM, characterized by post-glacial warming and continental ice-sheet retreat.

3.1. mtDNA Diversity

Arctic fox mitochondrial diversity was estimated using the complete dataset, comprising 80 complete mitogenomes (16,414 bp), representing 33 haplotypes, as well as geographic subsets thereof (Table 2). For the entire dataset, haplotype diversity is 0.887 ± 0.029, and nucleotide diversity is 0.0034 ± 0.0001. Among populations, Russia had the highest diversity (Hd = 1.00 ± 0.03, Pi = 0.005 ± 0.001), while samples from Fennoscandia had the lowest observed values (Hd = 0.65 ± 0.06, Pi = 0.0020 ± 0.0002). The average number of nucleotide differences between sequences within a population ranged from 85.33 (East Russia) to 17.33 (Norway), revealing a non-uniform distribution of mitochondrial diversity across the species’ range. However, we note that our results will be impacted by unequal sampling within regions.

3.2. mtDNA Diversity in Iceland

Among the 22 Icelandic samples, we observed 13 haplotypes (Hd = 0.939 ± 0.029; Pi = 0.00260 ± 0.00031; Table 2), none of which are shared with samples from elsewhere (Figure 2 and Figure 3), and that belong to three haplogroups (Hg.3, Hg.4, Hg.6; Figure 4). The two Icelandic haplotypes of Hg.3 are not closely related, differing by 38 substitutions (Figure 2). Three regions of Iceland (north-west, south-west, and north-east) harboured mitogenomes belonging to only a single haplogroup (Figure 4). The north-west coastal population has four haplotypes from Hg.6 (six samples, Hd = 0.80 ± 0.172, Pi = 8 × 10^−5^ ± 3 × 10^−5^), the south-west coastal population has one haplotype from Hg.4 (three samples, Hd = 0, Pi = 0), and the north-east inland population has one haplotype from Hg.3 (four samples, Hd = 0, Pi = 0). The southern inland population has eight haplotypes, distributed among Hg.3, Hg.4 and Hg.6 (nine samples, Hd = 0.97 ± 0.064, Pi = 0.0024 ± 0.00069).

When data are considered by ecotype, inland foxes showed slightly higher diversity (Hd = 0.910 ± 0.068; Pi = 0.00276 ± 0.00035) compared with coastal samples (Hd = 0.833 ± 0.098; Pi = 0.00218 ± 0.00054) (Table 2).

4. Discussion

The enhanced phylogenetic and temporal resolution provided by complete mitogenomes has improved inferences about the role of palaeoclimatic oscillations in matrilineal divergence and colonization routes in other circumpolar mammals. In brown bears, mitogenome data resolved five distinct clades not apparent with shorter mtDNA sequences and tied their divergence to post-glacial migration [28]. In reindeer, mitogenome data helped to resolve competing High Arctic Island colonization scenarios and linked this to early-Holocene range expansion [61]. In grey wolves, mitogenome data tied palaeoclimatic oscillations to the colonization of the Japanese islands and to intercontinental exchange between Eurasia and North America [29].

Similarly, increasing analyses of mitochondrial diversity from short fragments to whole mitochondrial genomes improved our ability to resolve matrilines in the Arctic fox. Our results refine the previous two-clade framework based on partial mtDNA sequences [26,27], increasing the number of matrilines found in the Holarctic to seven (Hg.1–Hg.7). These newly resolved haplogroups reveal a history of diversification shaped by alternating phases of glaciation and the expansion of open, shrub-dominated steppe-tundra across MIS 6-3 [59,62].

Three haplogroups (Hg.1, Hg.3, and Hg.5) are present in eastern Russia, two of which are not found elsewhere (Hg.1, Hg.5), which includes the most divergent haplogroup (Hg.1). These findings are consistent with a proposed refugial area in eastern Siberia during glacial maxima in the late Middle Pleistocene and Late Pleistocene [7,63], and indicate that two (extant) matrilines persisted through the Eemian interglacial in MIS 5e (~125 kya). Previous evidence for a matriline predating and surviving the Eemian interglacial comes from archival samples (partial mitogenomes [27]). Because our dataset included only a few central and eastern Russian samples, and notably only one Nearctic sample (from Canada), the full geographic distribution of Hg.1 and Hg.5 remains uncertain. Additionally, the existence of Nearctic matrilines that persisted through the Eemian cannot be excluded.

The two previously identified clades (designated A and B, or B and A, in Larsson et al. 2019; Panitsina et al. 2023 [26,27]) correspond to Hg.2 and Hg.3 in our study. Consistent with those studies, Hg.2 is observed in extant Arctic fox populations in Fennoscandia and central Russia, while Hg.3 is observed in extant populations in Fennoscandia, central and eastern Russia, Canada, and now also Iceland. It should be noted that our results indicate these lineages diverged ~123 kya (95% HPD = 69–184 kya), following the Eemian interglacial (in MIS 5e), when global temperatures started to decline. This differs from previous estimates of divergence at 65 kya [26] and 70 kya [27], in MIS 4, when Eurasian ice sheets extended almost as far south as during the LGM [57], and suggests a different demographic process for the divergence: post-interglacial range expansion vs. divergence in glacial refugia. A post-Eemian divergence is also observed in a previous partial-mitogenome phylogeny, involving a 30 ky old sample (PAF07) from eastern Russia [27].

A major radiation ~100–90 kya following the interstadial in MIS 5c gave rise to Hg.3, Hg.4, Hg.5, and the ancestral lineage leading to Hg.6 and Hg.7. It has been suggested that Arctic foxes extended their range in this post-interstadial period [27], likely during the expansion of steppe-tundra following the interstadial warm phase [59]. Subsequent branching of Hg.6 and Hg.7 during MIS 4, in contrast, likely reflects divergence in isolation during extensive glaciation when ice sheets and fragmented tundra would have restricted population connectivity. Later, additional divergences in Hg.2 and Hg.3 coincide with cooling following interstadial warming at the beginning of MIS 3. The youngest divergences (tip-level) occur after the LGM, likely reflecting post-glacial expansion. Together, these patterns suggest that isolation during glacial periods (e.g., MIS 4) as well as periods of expansion of suitable habitat for the species, either following interglacial/interstadial warming (e.g., MIS 5b, MIS 3) or post-glacial cooling (e.g., post-LGM), have repeatedly structured matrilineal diversity in Arctic foxes.

The timing of inferred divergence events here should be viewed as provisional. Molecular clock calibrations can vary considerably depending on dataset composition and analytical framework, and such differences can shift node ages and, in the case of the Arctic fox, alter their inferred association with interglacial/interstadial or glacial phases. Additional work, ideally involving broad sampling across the Nearctic and eastern Palearctic, will hopefully improve our understanding of the drivers of lineage formation in Arctic foxes. This would not only clarify the historical role of climatic oscillations in shaping Arctic fox diversity but also provide context for anticipating their responses to ongoing Arctic warming. Additionally, the absence of haplogroup representatives in specific regions (e.g., Palearctic) does not indicate geographic exclusivity and does not rule out the possibility that these lineages were formerly present in those regions but subsequently lost due to genetic drift or local extinction.

4.1. Colonization of Iceland

Considering that only a single Nearctic mitogenome is available, our phylogenetic analyses cannot resolve the geographic source of Iceland’s Arctic foxes, although current mitochondrial diversity in Iceland indicates the arrival of several unrelated females. It is uncertain whether the latter arrived contemporaneously or during successive dispersal episodes. A Nearctic origin of Icelandic Arctic foxes is a plausible scenario given available sampling, but alternative explanations involving unsampled Eurasian contributions cannot be excluded, nor can the possibility that lineages formerly present in the Palearctic contributed to the colonization of Iceland.

Two haplogroups (Hg.4 and Hg.6) are currently unique to Iceland, and their geographic origins cannot be determined with confidence. In contrast, the third haplogroup (Hg.3) is present in the Nearctic (Canada) and is also widespread across the Palearctic, indicating that both Nearctic and Palearctic contributions to Iceland are possible. A previous study using short mtDNA fragments identified links between Iceland and Nearctic regions through shared or closely related haplotypes [25]. Broader sampling across Greenland and Canada is necessary to evaluate whether Icelandic haplotypes represent unsampled Nearctic diversity, and to clarify the distribution of Hg.3 in the Nearctic.

While we identified three haplogroups in Iceland, the divergence among Icelandic Hg.3 haplotypes indicates that at least four unrelated females are required to account for haplotypic diversity on the island. Hg.3 and Hg.4 are represented by few haplotypes (two and three, respectively), but Hg.6 had sufficient internal diversity (eight haplotypes in twelve samples) to estimate the coalescence time for this haplogroup, giving a TMRCA of ~5600 years (95% HPD = 1.4–11 kya). This age suggests a mid-Holocene colonization following the Holocene Thermal Maximum [64,65], after which seasonal sea-ice and drift-ice export from Greenland resumed [66,67], potentially enabling dispersal to Iceland. It should be noted that Hg.6 is not star-like; multiple related haplotypes may have arrived on Iceland together, which would reduce the colonization age relative to our TMRCA estimate, as not all variation would have arisen in situ in such a scenario. Alternatively, colonization may have occurred earlier, during the early Holocene, while tundra habitat was present on Iceland [68], and reference therein), and before postglacial warming reduced connectivity between Greenland and Iceland. This earlier arrival would require a contraction of effective population size to produce the limited diversity observed, perhaps reflecting a bottleneck associated with forest expansion in Iceland until 6 kya [68].

4.2. Mitochondrial Diversity in Iceland

Haplogroups are not distributed equally in Iceland (Figure 4). Three regions are characterized by only one haplogroup each, in two cases by a single haplotype. Given the small per-site sample sizes and limited spatial coverage, additional haplotypes are likely present, and this pattern probably reflects incomplete detection rather than true fixation. The south-west coastal region has one haplotype from Hg.4 (3 samples), and the eastern inland region has one haplotype from Hg.3 (4 samples). The north-western coastal region has four haplotypes from Hg.6 (6 samples). The southern inland region is the most diverse, with representatives of all three haplogroups. Mitochondrial diversity is thus not structured by ecotype. Our results suggest that, at least historically, there has been gene flow between the western coastal region and the southern inland region. Results from nuclear microsatellite data [21] suggest that this has not been recent with regard to the north-western region. While the two inland regions have haplotypes belonging to Hg.3, they harbour different haplotypes that are quite distinct (38 substitutions); from this, we cannot surmise recent or ongoing (female) dispersal among these inland locations.

Why some regions harbour only a single haplogroup is yet unknown. Was there a ‘cosmopolitan’ (i.e., diverse) historical population, followed by drift in small, semi-isolated populations (e.g., during the 1970s bottleneck, or earlier)? Or does the diverse inland haplogroup composition reflect historical admixture after the initial coastal establishment of Arctic foxes in Iceland? To clarify the relative contributions of colonization history, drift, and isolation to Icelandic diversity, broader sampling is needed. Ideally, this should incorporate multiple marker types, as mitochondrial and nuclear markers provide complementary temporal perspectives, with mtDNA reflecting female dispersal and capturing deeper, colonization-related processes, while nuclear markers reflect more recent population structure.

5. Conclusions

Our findings show that Arctic foxes have repeatedly responded to climate oscillations, alternating between isolation and expansion as the circumpolar environment shifted. However, sparse sampling in key areas limits our ability to fully resolve these dynamics or to identify the geographic source of Iceland’s Arctic foxes. Integrating genomic and palaeoenvironmental data across the species’ range (especially from the Nearctic) will be crucial to understand how environmental fluctuations have shaped Arctic fox diversity, predicting how this species may respond to the current climate crisis [69].