Abstract

The circumpolar arctic fox Alopex lagopus thrives in cold climates and has a high migration rate involving long-distance movements. Thus, it differs from many temperate taxa that were subjected to cyclical restriction in glacial refugia during the Ice Ages. We investigated population history and genetic structure through mitochondrial control region variation in 191 arctic foxes from throughout the arctic. Several haplotypes had a Holarctic distribution and no phylogeographical structure was found. Furthermore, there was no difference in haplotype diversity between populations inhabiting previously glaciated and unglaciated regions. This suggests current gene flow among the studied populations, with the exception of those in Iceland, which is surrounded by year-round open water. Arctic foxes have often been separated into two ecotypes: ‘lemming’ and ‘coastal’. An analysis of molecular variance suggested particularly high gene flow among populations of the ‘lemming’ ecotype. This could be explained by their higher migration rate and reduced fitness in migrants between ecotypes. A mismatch analysis indicated a sudden expansion in population size around 118 000 BP, which coincides with the last interglacial. We propose that glacial cycles affected the arctic fox in a way opposite to their effect on temperate species, with interglacials leading to short-term isolation in northern refugia.

INTRODUCTION

The Quaternary cold periods are considered to have had a strong influence on the geographical distribution and genetic variation of organisms worldwide. In continental Eurasia and North America, repeated glaciations caused multiple periods of isolation in southern refugia and resulted in increased intraspecific genetic divergence (Taberlet et al., 1998; Hewitt, 2001). Several mammal species display phylogeographical patterns predicted by the expansion/contraction model with, for example, a high divergence between phylogroups from different refugia and genetic signatures of late Pleistocene expansions in population size (Hewitt, 1996). However, in highly mobile species gene flow during interglacials could lead to an admixture of genotypes from different refugia (Cruzan & Templeton, 2000). Furthermore, the impact of glaciation would have been different in species that were well adapted to cold climates compared with temperate species (Hewitt, 2001). Arctic species will not have been in southern temperate refugia and should thus not display the expansions/contractions associated with them. Arctic species may, however, have gone through range changes and they could have had different glacial and/or interglacial refugia.

The arctic fox Alopex lagopus is well adapted to arctic conditions (Fuglei & Øritsland, 1999) and in winter fur tolerates ambient temperatures below −40 °C without having to increase its metabolic rate significantly to keep a constant body temperature (Scholander et al., 1950). Its diet is composed of a variety of vertebrates (Audet, Robbins & Larivière, 2002), but two ecotypes are generally recognized: ‘lemming foxes’ that feed mainly on lemmings (Lemmus spp. and Dicrostonyx spp.) and ‘coastal foxes’ that feed mainly on eggs, birds and carrion from the marine system (Braestrup, 1941). Lemming foxes are found in continental Eurasia, North America, the Canadian archipelago and east Greenland, whereas coastal foxes are found in Iceland, Svalbard and south, west and north-west Greenland (Tannerfeldt & Angerbjörn, 1998). The difference between a highly fluctuating food source (lemming) and one that is more stable (coastal) has led to a number of different life-history strategies, where lemming foxes undergo an enormous reproductive output during lemming peaks compared with coastal foxes (Tannerfeldt & Angerbjörn, 1998). Furthermore, there are significant differences in migration patterns between the two ecotypes, with lemming foxes migrating further than coastal foxes (Angerbjörn, Hersteinsson & Tannerfeldt, 2004a).

Several studies suggest a high migration rate in A. lagopus, and that they are capable of long (> 1000 km) movements over the polar pack ice (e.g. Eberhardt & Hansson, 1978). Several subspecies of A. lagopus have been proposed, for example A. l. fuliginosus (Iceland), A. l. groenlandicus (Greenland), A. l. spitzbergenensis (Svalbard) and A. l. ungava (Canada) (Audet et al., 2002). Frafjord (1993) found some latitudinal differences in morphology between populations on a circumpolar scale, but pointed out that more information was needed on the genetic differentiation among A. lagopus populations.

In this study, we analysed mitochondrial DNA (mtDNA) variation in A. lagopus on a circumpolar scale to investigate the genetic structure and population history of the species. Concerning the population history, we did not expect to find the patterns of a rapid postglacial increase in population size which have been observed in more temperate species, since the wide distribution of A. lagopus during the last Ice Age (Kurtén, 1968; Kurtén & Anderson, 1980) suggests that A. lagopus were at least as abundant during this period as they are today. Instead, it is more probable that the warm interglacials have had a negative effect on the abundance of A. lagopus. We did, however, expect to see phylogeographical patterns from a postglacial range expansion in A. lagopus since they must have colonized formerly glaciated areas at the end of the last Ice Age. Past fragmentation events may be inferred from genetic distance among haplotypes, and their spatial distribution provides information on current gene flow among populations (Avise et al., 1987). Based on the high migration rate and long-distance movements observed in A. lagopus, we hypothesized that there is gene flow between most sampled populations that are connected via land or the polar sea ice (i.e. all populations except that in Iceland). We therefore predicted little phylogeographical structure, low ΦST values between all populations except that in Iceland and that populations in previously glaciated regions should have a haplotype diversity similar to those in continuously unglaciated regions. We also examined the long-term effective female population size and compared this with current estimates of the worldwide population size.

MATERIAL AND METHODS

DNA samples were collected from 191 A. lagopus from 13 regions throughout the arctic. The regions sampled were Svalbard (SVA), Iceland (ICE), east Greenland (EG), south Greenland (SG), west Greenland (WG), north-west Greenland (NWG), Churchill Manitoba (CHU), Cambridge Bay (CMB), Bathurst Island (BAT), Banks Island (BAN), Alaska (ALA), Siberia (SIB) and Fennoscandia (FEN) (Fig. 1). For statistical analyses concerned with geographical distances, we divided Siberia into east and west Siberia (two samples from Taimyr, which is halfway between east and west Siberia, were excluded from these analyses along with one sample for which it was not clear whether it was from east or west Siberia). Cambridge Bay, Banks Island and Bathurst Island were in some instances pooled into Canadian Archipelago (CA) in order to increase statistical power (there was no significant genetic differentiation among the regions within each pooling). Thirty-two of the samples were from the previous study by Dalén et al. (2002). The samples from Greenland were those previously used for microsatellite analysis by Meinke, Kapel & Arctander (2001). Tissue samples from Alaska were obtained from the University of Alaska Museum (UAM AF371–AF377, AF379, AF4012–AF4014, AF4039, AF21094).

Figure 1

Sample sites and number of samples (indicated within each circle) from each location. Populations clockwise from Greenwich Mean Time are: Iceland (ICE), east Greenland (EG), south Greenland (SG), west Greenland (WG), north-west Greenland (NWG), Churchill (CHU), Bathurst Island (BAT), Cambridge Bay (CMB), Banks Island (BAN), Alaska (ALA), Siberia (SIB), Fennoscandia (FEN) and Svalbard (SVA). The light grey area inside the dashed line illustrates the extent of polar sea ice in January (data from EOSDIS NSIDC Distributed Active Archive Center, http://nsidc.org/data/index.htm).

Figure 1

Sample sites and number of samples (indicated within each circle) from each location. Populations clockwise from Greenwich Mean Time are: Iceland (ICE), east Greenland (EG), south Greenland (SG), west Greenland (WG), north-west Greenland (NWG), Churchill (CHU), Bathurst Island (BAT), Cambridge Bay (CMB), Banks Island (BAN), Alaska (ALA), Siberia (SIB), Fennoscandia (FEN) and Svalbard (SVA). The light grey area inside the dashed line illustrates the extent of polar sea ice in January (data from EOSDIS NSIDC Distributed Active Archive Center, http://nsidc.org/data/index.htm).

Whole genomic DNA was extracted using Qiagen's Dneasy tissue kit (Qiagen). Faecal DNA (N = 4) was extracted from c. 200 mg dried faecal matter using the Qiaamp DNA stool mini kit (Qiagen). An approximately 320-bp fragment of the mitochondrial control region was amplified as previously described in Dalén et al. (2002). Sequencing of both the heavy and light strands was carried out using a CEQ 2000 automated sequencer (Beckman Coulter) following the manufacturer's instructions.

Sequences were aligned in BioEdit version 5.0.9 (Hall, 1999), checked by eye and assigned to haplotypes, which were named after their origin (see Fig. 2). We used the program ModelTest (Posada & Crandall, 1998) to evaluate which model of nucleotide substitution gave the best fit to the data. Sequence variability and population pairwise comparisons were computed with the software ARLEQUIN version 2.0 (Schneider, Roessli & Excoffier, 2000). Of the nucleotide substitution models supported in Arlequin, the Tamura & Nei (1993) model gave the lowest log likelihood score (with a gamma parameter of 0.7), and this was subsequently used in further analyses. Sequence variability was estimated as haplotype diversity (H), nucleotide diversity (πNei, 1987) and the mean number of pairwise differences (Tajima, 1993). Historic demographic expansions of population size were investigated through a mismatch analysis where the distribution of pairwise differences was compared with the expected distribution under a model of sudden expansion (Rogers & Harpending, 1992; Schneider & Excoffier, 1999). The estimated time of sudden expansion can be calculated from the equation τ = 2µt (Rogers, 1995), where µ is the mutation rate for the sequence and t is the time since expansion (confidence intervals for τ were obtained from 2000 bootstrap replicates). We also performed Fu's test of selective neutrality with 10 000 bootstrap replicates, where a significant negative Fs value (at P < 0.02) would indicate an expansion in population size (Fu, 1997). The long-term effective female population size (Nf) was approximated using the equation: Nf = 106(p/s)/g where p is the nucleotide diversity, s is the rate of sequence divergence and g is the generation time in years (Wilson et al., 1985). For the above calculations, we assumed a rate of sequence divergence of 14.2% Myr−1 (µ = 2.073 × 10−5), a rate that was recently estimated for wolves and coyotes (Savolainen et al., 2002), and a generation time of 2 years.

Figure 2

The minimum spanning network. Haplotypes are named after geographical origin: Holarctic (H), Nearctic (N), Canada (C), Siberia (S), Greenland (G) and Iceland (I). Each branch represents one mutational step; missing haplotypes are represented by a dot. Equally parsimonious branches are shown with dashed lines. The shape of the haplotypes illustrates the second nesting level in the nested clade analysis. Haplotype G3 was not nested until the third nesting level.

Figure 2

The minimum spanning network. Haplotypes are named after geographical origin: Holarctic (H), Nearctic (N), Canada (C), Siberia (S), Greenland (G) and Iceland (I). Each branch represents one mutational step; missing haplotypes are represented by a dot. Equally parsimonious branches are shown with dashed lines. The shape of the haplotypes illustrates the second nesting level in the nested clade analysis. Haplotype G3 was not nested until the third nesting level.

Using a function implemented in ARLEQUIN, we constructed a minimum spanning network based on pairwise differences among haplotypes (including indels). This network was subsequently used in a nested clade analysis (NCA) in an attempt to discriminate between phylogeographical patterns caused by the current restricted gene flow and patterns caused by historical events (Templeton, 1998). Nesting of the minimum spanning network followed the basic rules by Templeton, Boerwinkle & Sing (1987). Nesting of ambiguities and intermediate haplotypes was carried out according to Templeton & Sing (1993) and Crandall (1996). Geographical distances between regions were obtained using the distance calculator at http://www.wcrl.ars.usda.gov/cec/java/lat-long.htm (2003-01-28, Byers, 1997). The null hypothesis of no geographical associations of clades was tested and computation of clade distances and nested clade distances were carried out using the program GeoDis (Posada, Crandall & Templeton, 2000) with 10 000 permutations. Interpretation of the results obtained in the NCA was obtained using the inference key in GeoDis.

We employed an exact non-parametric procedure (1 000 000 steps in the Markov chain and 50 000 dememorization steps) to test for differentiation between pairs of populations (Raymond & Rousset, 1995). In order to investigate geographical structuring of genetic variation, we used an analysis of molecular variance (AMOVA) with 10 000 permutations (Excoffier, Smouse & Quattro, 1992). We performed six AMOVAs with different hierarchical groupings: [Palaearctic vs. Nearctic], [Palaearctic vs. Nearctic vs. Atlantic islands], [mainland vs. islands], [above 68°N vs. below 68°N], [lemming fox populations vs. coastal fox populations] and [lemming fox populations vs. each coastal fox population]. We then assumed that the most probable geographical structure was represented by the groupings that maximized values of ΦCT (Vila et al., 1999), which is a measure of the proportion of genetic variation among groupings of populations. Population pairwise ΦST values (a measure analogous to FST) were generated and tested for significance through 10 000 permutations (Schneider et al., 2000). The resulting matrix of ΦST values between the different populations was visualized with a UPGMA tree constructed in PAUP (Swofford, 1998). To investigate the effect of postglacial gene flow, we compared H for populations in formerly glaciated areas with those inhabiting regions not glaciated during the last Ice Age. This was done with a one-way ANOVA, as implemented in the software STATISTICA (StatSoft Inc., 1999). For this analysis we excluded samples from Bathurst Island due to low sample size. A Mantel test with 10 000 replicates (Smouse, Long & Sokal, 1986) was used to test if there was a correlation between genetic and geographical distances among populations.

RESULTS

We sequenced 292 bp of the control region for each of the 191 individuals. The sequenced region contained 21 variable sites, which defined 29 different haplotypes (Table 1). All the observed variation was in the form of single base-pair substitutions or indels, except for haplotype S3, in which a 16-bp deletion was observed (since this region was present in all other haplotypes, as well as in kit and swift foxes (Vulpes macrotis and V. velox), it was presumably a deletion). This deletion was confirmed by a second amplification and sequence analysis using two additional primers, H1F (5′-GCCATCAACTCCCAAAGCT-3′) and P1R (5′-GAGGCATGGTGATAAATCC-3′). The whole deletion was treated with the same weight as substitutions and indels in further statistical analyses. The mean number of pairwise differences between all samples was 2.65 (SD, 1.42), and π in the total sample was 0.009 (SD, 0.005). Fu's test of selective neutrality gave a significantly large negative Fs value (Fs = −8.15, P = 0.014).

Table 1

Geographical distribution and GenBank accession numbers for Alopex lagopus haplotypes

Haplotype GenBank ♯ Geographical region
 
FEN SIB ICE BAT CHU CMB WG EG NWG SG ALA SVA BAN 
H1 AY321121 12    
H2 AY321125     
H3 AY321120           
H4 AY321124           
H5 AY321127          
H6 AY321128            
H7 AY321129           
H8 AY321132          
H9 AY321134           
N1 AY321136           
N2 AY321138           
N3 AY321140            
S1 AY321123             
S2 AY321133             
S3 AY321122             
S4 AY321126             
I1 AY321131             
I2 AY321130   14           
G1 AY321135             
G2 AY321137             
G3 AY321139             
G4 AY321141             
G5 AY321142             
C1 AY321143            
C2 AY321144             
C3 AY321145           
C4 AY321146             
C5 AY321147             
C6 AY321148             
Total  22 25 23 20 15 11 10 10 13 20 10 
 0.70 0.76 0.61 – 0.76 0.84 0.58 0.76 0.47 0.84 0.91 0.76 0.89 
Haplotype GenBank ♯ Geographical region
 
FEN SIB ICE BAT CHU CMB WG EG NWG SG ALA SVA BAN 
H1 AY321121 12    
H2 AY321125     
H3 AY321120           
H4 AY321124           
H5 AY321127          
H6 AY321128            
H7 AY321129           
H8 AY321132          
H9 AY321134           
N1 AY321136           
N2 AY321138           
N3 AY321140            
S1 AY321123             
S2 AY321133             
S3 AY321122             
S4 AY321126             
I1 AY321131             
I2 AY321130   14           
G1 AY321135             
G2 AY321137             
G3 AY321139             
G4 AY321141             
G5 AY321142             
C1 AY321143            
C2 AY321144             
C3 AY321145           
C4 AY321146             
C5 AY321147             
C6 AY321148             
Total  22 25 23 20 15 11 10 10 13 20 10 
 0.70 0.76 0.61 – 0.76 0.84 0.58 0.76 0.47 0.84 0.91 0.76 0.89 

The number of samples and haplotype diversities (H) are indicated for each geographical region (abbreviations as in Fig. 1).

The distribution of pairwise differences between all individuals did not deviate from the expected distribution under a model of sudden expansion (P = 0.45). The extent of divergence was measured as τ  = 4.889 (95% CI, 1.674–9.298), giving an estimated time of expansion at 118 000 BP (95% CI, 40 000–224 000). The observed nucleotide diversity suggested a long-term female effective population size of 32 000 (± 17 000) individuals in the sampling area (i.e. the world population).

The minimum spanning network of the different haplotypes revealed no major branching events (Fig. 2). Two haplotypes, H1 and H2, were observed in 42% of all individuals. These two haplotypes, together with several less common haplotypes, had a widespread geographical distribution. The remaining haplotypes were generally site-specific and occurred in low frequencies (Table 1). Haplotypes specific to certain geographical regions did not form monophyletic groups but instead appeared to be randomly distributed in the network (Fig. 2). The NCA did, however, indicate a significant geographical association for a majority of the nested clades (data not shown). We inferred that the overall phylogeographical pattern in the NCA was caused by recurrent but restricted gene flow. This pattern was dominant at the second and third (total network) nesting levels. At the first nesting level the pattern was more complicated, with indications of past fragmentations, range expansions and restricted gene flow (see Appendix for a complete listing of NCA results).

Figure 3

Population tree based on ΦST values, illustrating the most probable geographical structure in the analysis of molecular variance (AMOVA). The results suggest that there is high gene flow between populations belonging to the lemming ecotype, whereas gene flow seems to be lower between populations of the coastal ecotype as well as between the two ecotypes.

Figure 3

Population tree based on ΦST values, illustrating the most probable geographical structure in the analysis of molecular variance (AMOVA). The results suggest that there is high gene flow between populations belonging to the lemming ecotype, whereas gene flow seems to be lower between populations of the coastal ecotype as well as between the two ecotypes.

There was no significant correlation between genetic and geographical distances among populations (Mantel test: r = −0.19, P = 0.90). The exact test of population differentiation indicated that most populations were differentiated although there were exceptions, especially within North America (Table 2). The most probable geographical grouping of populations in the AMOVA was when lemming fox populations were grouped against each of the coastal fox populations (P < 0.002), where 25.4% of the variation was observed among groups (ΦCT values for other groupings were all below 3%). The total proportion of variation among all populations (ΦST) was 30%, and the proportion of variation among populations within groups (ΦSC) was 6.8%. Among the populations, the ΦST values were generally low with the exception of Iceland and to some extent west Greenland (Table 2). The H in the different populations varied between 0.47 and 0.91 (Table 1). There was no significant difference  in H between previously glaciated and non-glaciated regions (one way ANOVA, N = 13, F = 2.44, P = 0.15), where Banks Island, Alaska and East Siberia were considered as having been unglaciated during the latest Ice Age.

Table 2

Population differentiation test (P-values; above diagonal) and cross-wise ΦST values for each population (below diagonal)

 FEN SIB ICE CHU CA WG EG NWG SG ALA SVA 
FEN  0.030 <0.001 <0.001 <0.001 <0.001 0.001 <0.001 <0.001 <0.001 <0.001 
SIB 0.116  <0.001 0.021 0.010 <0.001 0.029 <0.001 0.002 0.029 0.001 
ICE 0.280 0.365  <0.001 <0.001 <0.001 <0.001 <0.001 <0.001 <0.001 <0.001 
CHU 0.195 0.025 0.467  0.557 <0.001 0.235 0.001 <0.001 0.079 0.010 
CA 0.080 0.007 0.357 0.008  0.001 0.343 0.002 0.003 0.193 0.024 
WG 0.390 0.417 0.537 0.555 0.419  0.001 <0.001 0.001 0.055 0.046 
EG 0.094 0.004 0.410 0.084 0.028 0.399  0.053 0.002 0.349 0.030 
NWG 0.361 0.197 0.564 0.397 0.286 0.505 0.144  <0.001 0.004 <0.001 
SG 0.135 0.183 0.407 0.282 0.130 0.256 0.219 0.466  0.004 <0.001 
ALA 0.043 0.074 0.204 0.179 0.052 0.210 0.076 0.272 0.063  0.026 
SVA 0.187 0.110 0.439 0.201 0.101 0.155 0.088 0.216 0.090 0.063  
 FEN SIB ICE CHU CA WG EG NWG SG ALA SVA 
FEN  0.030 <0.001 <0.001 <0.001 <0.001 0.001 <0.001 <0.001 <0.001 <0.001 
SIB 0.116  <0.001 0.021 0.010 <0.001 0.029 <0.001 0.002 0.029 0.001 
ICE 0.280 0.365  <0.001 <0.001 <0.001 <0.001 <0.001 <0.001 <0.001 <0.001 
CHU 0.195 0.025 0.467  0.557 <0.001 0.235 0.001 <0.001 0.079 0.010 
CA 0.080 0.007 0.357 0.008  0.001 0.343 0.002 0.003 0.193 0.024 
WG 0.390 0.417 0.537 0.555 0.419  0.001 <0.001 0.001 0.055 0.046 
EG 0.094 0.004 0.410 0.084 0.028 0.399  0.053 0.002 0.349 0.030 
NWG 0.361 0.197 0.564 0.397 0.286 0.505 0.144  <0.001 0.004 <0.001 
SG 0.135 0.183 0.407 0.282 0.130 0.256 0.219 0.466  0.004 <0.001 
ALA 0.043 0.074 0.204 0.179 0.052 0.210 0.076 0.272 0.063  0.026 
SVA 0.187 0.110 0.439 0.201 0.101 0.155 0.088 0.216 0.090 0.063  

ΦST values not significantly different from zero. Population abbreviations as in Fig. 1.

DISCUSSION

Population history

The earliest historic event that can be inferred from the mitochondrial DNA variation is that of a sudden expansion in population size. This was presumably preceded by a population bottleneck. The occurrence of a historic bottleneck and subsequent expansion is further supported by the significantly negative Fs value (Fu, 1997) and by the low nucleotide diversity (0.009) in A. lagopus. The nucleotide diversity in the control region was considerably lower than that in other mammals, for example wolves (Canis lupus; π  = 0.026; Vila et al., 1999), coyotes (Canis latrans; π = 0.046; Vila et al., 1999) and moose (Alces alces; π = 0.025; Hundertmark et al., 2002). The time of the expansion, as suggested by the mismatch analysis, was estimated at approximately 118 000 BP. Bearing in mind the large confidence interval (40 000–224 000 BP), this estimate coincides with the last interglacial which ended in 117 000 BP (Kukla et al., 2002). A similar expansion in connection with the last interglacial was recently observed in a study on reindeer (Flagstad & Røed, 2003). Considering that the last interglacial was approximately 5 °C warmer than temperatures are at present (Funder et al., 1998), it is probable that A. lagopus (along with other arctic organisms, such as reindeer) was adversely affected during this period. This may, for example, have been through indirect effects, such as a northern expansion of the red fox Vulpes vulpes, as it has been proposed that the southern distribution of A. lagopus is limited by V. vulpes (Hersteinsson & Macdonald, 1992). The presence of forest remains from previous interglacials in northern Siberia (Sher, 1991) suggests a suitable habitat for V. vulpes. A. lagopus may therefore have been extinct in continental Eurasia and North America during the last interglacial, persisting only in high-latitude islands, and then expanding south as temperatures started to fall some 117 000 years ago. This hypothesis predicts a high current genetic diversity in high-latitude islands that were not glaciated during the ensuing Ice Age, since these are the only areas that would have been continuously inhabited by A. lagopus for at least 130 000 years (the low sample sizes from these islands in our study did not allow us to test this hypothesis). It can also be expected that any sequences recovered from fossil remains less than 100 000 years old would fall within the scope of the mismatch distribution.

During the Ice Age that followed, A. lagopus was widely distributed in Eurasia and Beringia (Kurtén, 1968; Kurtén & Anderson, 1980). The structure of the minimum spanning network, without distinct phylogroups, indicates a lack of significant geographical barriers during this period (Fig. 2). This is further supported by the lack of evidence of past fragmentation at the higher nesting levels in the NCA (Appendix). At the end of the Ice Age it is probable that there was a range expansion into formerly glaciated areas such as Greenland, Iceland, Svalbard and Fennoscandia. Although the second nesting level in the NCA showed some support for a range expansion, this was not as strong as might have been expected. There could be several explanations for this; for example, that these areas were colonized by A. lagopus from local refugia (e.g. Frafjord & Hufthammer, 1994), as has been suggested for other arctic species (Fedorov & Stenseth, 2001, 2002). A high postglacial gene flow could also explain the weak support in the NCA as it may have erased phylogeographical patterns created by an initial range expansion. At the lowest nesting level, the NCA gave a rather ambiguous picture, possibly due to small sample sizes in the nested clades. There may also be problems with the interpretation of NCA results using the inference key (see Knowles & Maddison, 2002).

The female long-term effective population size was estimated at 32 000 individuals. Assuming a 1 : 1 sex ratio and that 40% of all female adult A. lagopus breed during their lifetime (Angerbjörn et al., 2004a), this would correspond to an approximate world population size of c. 160 000 (± 85 000) adults. This is lower than the census population size of 330 000–930 000 adults (Angerbjörn, Hersteinsson & Tannerfeldt, 2004b), but is within the margins of what can be expected for a species with a large variance in reproductive success (Creel, 1998; Bensch & Hasselquist, 1999). Thus, we did not find any indication of recent changes in the world-wide population size of A. lagopus as have been reported for other canids (e.g. Vila et al., 1999).

Current genetic structure

Although most populations seemed to be significantly differentiated from each other, several analyses suggested that currently there is restricted gene flow between the majority of the populations. There was no phylogeographical structure in the minimum spanning network, where presumed ancestral haplotypes were frequent and widespread and newly arisen haplotypes have not yet spread throughout the range of the species. Therefore, A. lagopus appears to be a species with intermediate gene flow and no long-term zoogeographical barriers (category V in Avise et al., 1987). A similar lack of phylogeographical structure has previously been observed in fish (e.g. Rocha-Olivares, Garber & Stuck, 2000), and to some extent wolves (Canis lupus; Vila et al., 1999). The predominantly low ΦST values among populations on such a large geographical scale, compared with ΦST values of 0.75 in kit foxes, 0.50 in swift foxes (Mercure et al., 1993), 0.46 in Mediterranean V. vulpes (Frati et al., 1998) and 0.69 in wolves (Vila et al., 1999), also indicate current gene flow. Yet low ΦST values and poor phylogeographical structuring of haplotypes could also be the result of a postglacial range expansion. However, the high ΦST values between Iceland and the other populations suggest that Iceland is particularly isolated, as would be expected under the hypothesis that there is current gene flow between all populations except that in Iceland. The observation of equal haplotype diversities in populations inhabiting formerly glaciated and unglaciated areas further supports the gene flow hypothesis, although it should be noted that colonization of a formerly glaciated region from several different refugia can also result in high haplotype diversity (Hewitt, 1996). Taken together, these results suggest that there is gene flow among most populations, which is in agreement with previous studies reporting that A. lagopus travels long distances (e.g. Eberhardt & Hanson, 1978) and illustrates the importance of the polar sea ice for terrestrial arctic mammals.

We could not, however, find a correlation between genetic and geographical distances, implying that there is no genetic isolation by distance between the populations. There could be a number of explanations for this, such as ice movements, geographical barriers or A. lagopus following polar bears (however, we could find no relationship between A. lagopus and polar bear genetic distances; Paetkau et al., 1999). A more likely explanation can be found in the relationship between the different populations and the geographical structuring of the genetic variation as suggested by the AMOVA. A. lagopus from east Greenland, Siberia, Churchill, the Canadian Archipelago, Fennoscandia and Alaska formed a group of populations that were genetically more closely related to each other than to any of the other populations. This former group consisted of populations with lemming foxes, whereas the latter populations were all of the coastal fox ecotype. As indicated by the AMOVA, only 6.8% of the genetic variation could be explained by differences among lemming fox populations whereas 25.4% of the variation could be explained by differences between the lemming fox group and each of the coastal fox populations. It therefore seems that gene flow is substantially higher between populations of lemming foxes than it is between the two ecotypes or between coastal fox populations. The ecological causes for such a pattern could be that lemming foxes have a higher frequency of long-distance migrations (Angerbjörn et al., 2004a), and that migrants from one type of habitat to the other have a lower fitness compared with resident A. lagopus. That lemming foxes should migrate longer and more often than do coastal foxes actually makes evolutionary sense owing to the large-scale spatial synchrony of lemming populations (Krebs et al., 2002), which may force foxes feeding on lemmings to migrate longer and more frequently than do foxes in coastal areas where food resources are more stable. The hypothesis that immigrant foxes from a different habitat should have lower fitness compared with residents was originally proposed by Vibe (1967) as an explanation for the stable difference in fur colour frequency between A. lagopus in north-west Greenland and Canada, despite an influx of white foxes after lemming peaks in Canada. It has been suggested that different reaction norms in litter size have evolved in fluctuating and stable A. lagopus populations (Tannerfeldt & Angerbjörn, 1998). The observed pattern might thus be explained if food resource predictability affects selection pressure on reproductive output, giving lemming foxes a disadvantage under stable coastal conditions, or by the higher competition for territories in coastal fox populations (Angerbjörn et al., 2004a).

These results agree well with what is known on the biology of A. lagopus, in particular the extraordinary migration patterns facilitated by the polar sea ice, and the difference in life-history strategies between lemming and coastal A. lagopus. The generally high gene flow suggested by this study, in particular among lemming fox populations, should also be taken into account with respect to the spread of arctic disease, such as rabies.

Management implications

Our samples covered more or less the total distribution of A. lagopus except for the populations on the isolated Bering and Mednyi Islands. We found no support for the existence of any subspecies within the sampled area. Furthermore, based on the distribution of mtDNA haplotypes, we were unable to identify any Evolutionary Significant Units. Iceland may, however, be considered a Management Unit based on its isolation, as indicated by the high ΦST values. However, Management Units should not be based solely on genetic data. Fennoscandia, for example, is regarded as a Management Unit based on ecological data. In a previous study by Dalén et al. (2002) it was suggested that there is gene flow from Siberia into Fennoscandia, since the haplotype diversity and number of haplotypes in Fennoscandia was higher than expected for a small isolated population. Two observations in this study support that conclusion. First, the ΦST values between Fennoscandia and Siberia (0.12) was not particularly high compared with the difference between other populations. Second, the two haplotypes that had previously been observed only in Fennoscandia were in this study also found in western Siberia, which is to be expected if the haplotypes in Fennoscandia are the result of current gene flow from Siberia.

On a global level, the results of this study suggest that the high temperatures during the last interglacial may have had a severe impact on A. lagopus as a species. Given the increases in temperature predicted from models on global warming and the negative effect of competition with the temperate V. vulpes (Chirkova, 1968; Tannerfeldt, Elmhagen & Angerbjörn, 2002), the range of A. lagopus will contract to the north. The local conservation problems for A. lagopus in Fennoscandia today may thus, in the near future, become a global issue.

ACKNOWLEDGEMENTS

We are grateful to Torsten Eriksson for help with phylogenetic analyses, and Anna Linderholm and Annica Olsson for help with primers. Olga Pavlova contributed with information on the extent of the polar sea ice. The University Museum of Alaska provided tissue samples from Alaska. We are also very grateful to all field personnel who helped to collect tissue samples, in particular Petteri Polojärvi and Risto Karvonen who provided samples from the Kola Peninsula, and Harald Solheim who trapped foxes in Svalbard. We are thankful to Øystein Flagstad for providing valuable comments on the manuscript. The Swedish Polar Research Secretariat organized the transpolar expeditions Tundra Ecology-94 and Tundra North-west 1999. All genetic analyses were financed by the Ebba & Sven Schwartz Foundation.

REFERENCES

Angerbjörn
A
,
Hersteinsson
P
,
Tannerfeldt
M
2004a
.
Consequences of resource predictability in the arctic fox – two life history strategies
. In:
Macdonald
DW
Sillero-Zubiri
C
eds.
The biology and conservation of wild canids
 .
Oxford
:
Oxford University Press
.
Angerbjörn
A
,
Hersteinsson
P
,
Tannerfeldt
M
2004b
.
Arctic fox
. In:
Macdonald
DW
Sillero-Zubiri
C
eds.
Canid action plan
 .
Gland
:
IUCN.
Audet
AM
,
Robbins
BC
,
Larivière
S
2002
.
Alopex lagopus
.
Mammalian Species
 
713
:
1
10
.
Avise
JC
,
Arnold
J
,
Ball
RM
,
Bermingham
E
,
Lamb
T
,
Neigel
JE
,
Reeb
C
,
Saunders
NC
1987
.
Intraspecific phylogeography – the mitochondrial-DNA bridge between population genetics and systematics
.
Annual Review of Ecology and Systematics
 
18
:
489
522
.
Bensch
S
,
Hasselquist
D
1999
.
Phylogeographic population structure of great reed warblers: an analysis of mtDNA control region sequences
.
Biological Journal of the Linnean Society
 
66
:
171
185
.
Braestrup
FW
1941
.
A study on the arctic fox in Greenland
.
Meddelelser om Grønland
 
13
(4) :
1
101
.
Byers
JA
1997
.
Surface distance between points of latitude and longitude
 . http://www.wcrl.ars.usda.gov/cec/java/lat-long.htm (2003-01-28).
Chirkova
AF
1968
.
A relationship between arctic fox and red fox in the far north
.
Problems of the North
 
11
:
129
131
.
Crandall
KA
1996
.
Multiple interspecies transmissions of human and simian T-cell leukemia/lymphoma virus type I sequences
.
Molecular Biology and Evolution
 
13
:
115
131
.
Creel
S
1998
.
Social organisation and effective population size in carnivores
. In:
Caro
T
ed.
Behavioral ecology and conservation biology
 .
NewYork
:
Oxford University Press.
Cruzan
MB
,
Templeton
AR
2000
.
Paleoecology and coalescence: phylogeographic analysis of hypotheses from the fossil record
.
Trends in Ecology and Evolution
 
15
:
491
496
.
Dalén
L
,
Götherström
A
,
Tannerfeldt
M
,
Angerbjörn
A
2002
.
Is the endangered Fennoscandian arctic fox (Alopex lagopus) population genetically isolated?
Biological Conservation
 
105
:
171
178
.
Eberhardt
L
,
Hansson
WC
1978
.
Long distance movements of arctic foxes tagged in northern Alaska
.
Canadian Field Naturalist
 
92
:
386
389
.
Excoffier
L
,
Smouse
P
,
Quattro
J
1992
.
Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data
.
Genetics
 
131
:
479
491
.
Fedorov
VB
,
Stenseth
NC
2001
.
Glacial survival of the Norwegian lemming (Lemmus lemmus) in Scandinavia: inference from mitochondrial DNA variation
.
Proceedings of the Royal Society of London B
 
268
:
809
814
.
Fedorov
VB
,
Stenseth
NC
2002
.
Multiple glacial refugia in the North American Arctic: inference from phylogeography of the collared lemming (Dicrostonyx groenlandicus)
.
Proceedings of the Royal Society of London B
 
269
:
2071
2077
.
Flagstad
Ø
,
Røed
KH
2003
.
Refugial origins of reindeer (Rangifer tarandus L.) inferred from mitochondrial DNA sequences
.
Evolution
 
57
:
658
670
.
Frafjord
K
1993
.
Circumpolar size variation in the skull of the arctic fox Alopex lagopus
.
Polar Biology
 
13
:
235
238
.
Frafjord
K
,
Hufthammer
AK
1994
.
Subfossil records of the arctic fox (Alopex lagopus) compared to its present distribution in Norway
.
Arctic
 
47
:
65
68
.
Frati
F
,
Hartl
GB
,
Lovari
S
,
Delibes
M
,
Markov
G
1998
.
Quaternary radiation and genetic structure of the red fox Vulpes vulpes in the Mediterranean Basin, as revealed by allozymes and mitochondrial DNA
.
Journal of Zoology
 
245
:
43
51
.
Fu
YX
1997
.
Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection
.
Genetics
 
147
:
915
925
.
Fuglei
E
,
Øritsland
NA
1999
.
Seasonal trends in body mass, food intake and resting metabolic rate, and induction of metabolic depression in arctic foxes (Alopex lagopus) at Svalbard
.
Journal of Comparative Physiology B
 
169
:
361
369
.
Funder
S
,
Hjort
C
,
Landvik
JY
,
Nam
SI
,
Reeh
N
,
Stein
R
1998
.
History of a stable ice margin in East Greenland during the middle and upper Pleistocene
.
Quaternary Science Reviews
 
17
:
77
123
.
Hall
TA
1999
.
BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT
.
Nucleic Acids Symposium
 
41
:
95
98
.
Hersteinsson
P
,
Macdonald
DW
1992
.
Interspecific competition and the geographical distribution of red and arctic foxes Vulpes vulpes and Alopex lagopus
.
Oikos
 
64
:
505
515
.
Hewitt
GM
1996
.
Some genetic consequences of ice ages, and their role in divergence and speciation
.
Biological Journal of the Linnean Society
 
58
:
247
276
.
Hewitt
G
2001
.
Speciation, hybrid zones and phylogeography – or seeing genes in space and time
.
Molecular Ecology
 
10
:
537
549
.
Hundertmark
KJ
,
Shields
GF
,
Udina
IG
,
Bowyer
RT
,
Danilkin
AA
,
Schwartz
CC
2002
.
Mitochondrial phylogeography of moose (Alces alces): late Pleistocene divergence and population expansion
.
Molecular Phylogenetics and Evolution
 
22
:
375
387
.
Knowles
LL
,
Maddison
WP
2002
.
Statistical phylogeography
.
Molecular Ecology
 
11
:
2623
2635
.
Krebs
CJ
,
Kenney
AJ
,
Gilbert
S
,
Danell
K
,
Angerbjörn
A
,
Erlinge
S
,
Bromley
RG
,
Shank
C
,
Carriere
S
2002
.
Synchrony in lemming and vole populations in the Canadian Arctic
.
Canadian Journal of Zoology
 
80
:
1323
1333
.
Kukla
GJ
,
Bender
ML
,
De Beaulieu
JL
,
Bond
G
,
Broecker
WS
,
Cleveringa
P
,
Gavin
JE
,
Herbert
TD
,
Imbrie
J
,
Jouzel
J
,
Keigwin
LD
,
Knudsen
KL
,
McManus
JF
,
Merkt
J
,
Muhs
DR
,
Muller
H
,
Poore
RZ
,
Porter
SC
,
Seret
G
,
Shackleton
NJ
,
Turner
C
,
Tzedakis
PC
,
Winograd
IJ
2002
.
Last interglacial climates
.
Quaternary Research
 
58
:
2
13
.
Kurtén
B
1968.
Pleistocene mammals of Europe
 .
London
:
Weidenfeld and Nicolson.
Kurtén
B
,
Anderson
E
1980.
Pleistocene mammals of North America
 .
NewYork
:
Columbia University Press
.
Meinke
PG
,
Kapel
CMO
,
Arctander
P
2001
.
Genetic differentiation of populations of Greenlandic Arctic Fox
.
Polar Research
 
20
:
75
83
.
Mercure
A
,
Ralls
K
,
Koepfli
KP
,
Wayne
RK
1993
.
Genetic subdivisions among small canids – mitochondrial-DNA differentiation of swift, kit, and arctic foxes
.
Evolution
 
47
:
1313
1328
.
Nei
M
1987.
Molecular evolutionary genetics
 .
NewYork
:
Columbia University Press
.
Paetkau
D
,
Amstrup
SC
,
Born
EW
,
Calvert
W
,
Derocher
AE
,
Garner
GW
,
Messier
,
F
,
Stirling
I
,
Taylor
MK
,
Wiig
O
,
Strobeck
C
1999
.
Genetic structure of the world's polar bear populations
.
Molecular Ecology
 
8
:
1571
1584
.
Posada
D
,
Crandall
KA
1998
.
Modeltest: testing the model of DNA substitution
.
Bioinformatics
 
14
:
817
818
.
Posada
D
,
Crandall
KA
,
Templeton
AR
2000
.
GeoDis: a program for the cladistic nested analysis of the geographical distribution of genetic haplotypes
.
Molecular Ecology
 
9
:
487
488
.
Raymond
M
,
Rousset
F
1995
.
An exact test for population differentiation
.
Evolution
 
49
:
1280
1283
.
Rocha-Olivares
A
,
Garber
NM
,
Stuck
KC
2000
.
High genetic diversity, large inter-oceanic divergence and historical demography of the striped mullet
.
Journal of Fish Biology
 
57
:
1134
1149
.
Rogers
AR
1995
.
Genetic evidence for a Pleistocene population explosion
.
Evolution
 
49
:
608
615
.
Rogers
AR
,
Harpending
H
1992
.
Population growth makes waves in the distribution of pairwise genetic differences
.
Molecular Biology and Evolution
 
9
:
552
569
.
Savolainen
P
,
Zhang
Y
,
Luo
J
,
Lundeberg
J
,
Leitner
T
2002
.
Genetic evidence for an East Asian origin of domestic dogs
.
Science
 
298
:
1610
1613
.
Schneider
S
,
Excoffier
L
1999
.
Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates vary among sites: application to human mitochondrial DNA
.
Genetics
 
152
:
1079
1089
.
Schneider
S
,
Roessli
D
,
Excoffier
L
2000.
Arlequin, Version 2.000. A software for population genetics data analysis
 .
Switzerland
:
Genetics and Biometry Laboratory, University of Geneva.
Scholander
PF
,
Hock
R
,
Walters
V
,
Johnson
F
,
Irving
L
1950
.
Heat regulation in some arctic and tropical mammals and birds
.
Biological Bulletin
 
99
:
237
258
.
Sher
AV
1991
.
Problems of the interglacial in Arctic Siberia
.
Quaternaty International
 
10–12
:
215
222
.
Smouse
PE
,
Long
JC
,
Sokal
RR
1986
.
Multiple regression and correlation extensions of the Mantel Test of matrix correspondence
.
Systematic Zoology
 
35
:
627
632
.
Swofford
DL
1998.
PAUP* phylogenetic analysis using parsimony (*and other methods).
 
Sunderland, MA
:
Sinauer Associates.
Taberlet
P
,
Fumagalli
L
,
Wust-Saucy
AG
,
Cosson
JF
1998
.
Comparative phylogeography and postglacial colonization routes in Europe
.
Molecular Ecology
 
7
:
453
464
.
Tajima
F
1993
.
Measurement of DNA polymorphism
. In:
Takhata
N
Clark
AG
eds.
Mechanisms of molecular evolution. Introduction to molecular paleopopulation biology.
  Tokyo, Japan and Sunderland, MA: Scientific
Societies Press, Sinauer Associates, Inc.
Tamura
K
,
Nei
M
1993
.
Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees
.
Molecular Biology and Evolution
 
10
:
526
.
Tannerfeldt
M
,
Angerbjörn
A
1998
.
Fluctuating resources and the evolution of litter size in the arctic fox
.
Oikos
 
83
:
545
559
.
Tannerfeldt
M
,
Elmhagen
B
,
Angerbjörn
A
2002
.
Exclusion by interference competition? The relationship between red and arctic foxes
.
Oecologia
 
132
:
213
220
.
Templeton
AR
1998
.
Nested clade analyses of phylogeographic data: testing hypotheses about gene flow and population history
.
Molecular Ecology
 
7
:
381
397
.
Templeton
AR
,
Boerwinkle
E
,
Sing
CF
1987
.
A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping. I. Basic theory and an analysis of alcohol dehydrogenase activity in Drosophila
.
Genetics
 
117
:
343
351
.
Templeton
AR
,
Sing
CF
1993
.
A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping. IV. Nested analyses with cladogram uncertainty and recombination
.
Genetics
 
134
:
659
669
.
Vibe
C
1967
.
Arctic animals in relation to climatic fluctuations. The arctic fox
.
Meddelelser om Grønland
 
170
:
101
150
.
Vila
C
,
Amorim
IR
,
Leonard
JA
,
Posada
D
,
Castroviejo
J
,
Petrucci-Fonseca
F
,
Crandall
KA
,
Ellegren
H
,
Wayne
RK
1999
.
Mitochondrial DNA phylogeography and population history of the grey wolf Canis lupus
.
Molecular Ecology
 
8
:
2089
2103
.
Wilson
AC
,
Cann
RI
,
Carr
SM
,
George
M
,
Gyllensten
U
,
Helm-Bychowski
KM
,
Higuchi
RG
,
Palumbi
SR
,
Prager
EM
,
Sage
RD
,
Stoneking
M
1985
.
Mitochondrial DNA and two perspectives on evolutionary genetics
.
Biological Journal of the Linnean Society
 
26
:
375
400
.

Appendix

Results of the nested clade analysis of geographical distances for control region haplotypes in Alopex lagopus

 
formula