Skip to main content

Advertisement

Fine-scale population genetic structure of arctic foxes (Vulpes lagopus) in the High Arctic

Abstract

Objective

The arctic fox (Vulpes lagopus) is a circumpolar species inhabiting all accessible Arctic tundra habitats. The species forms a panmictic population over areas connected by sea ice, but recently, kin clustering and population differentiation were detected even in regions where sea ice was present. The purpose of this study was to examine the genetic structure of a population in the High Arctic using a robust panel of highly polymorphic microsatellites.

Results

We analyzed the genotypes of 210 individuals from Bylot Island, Nunavut, Canada, using 15 microsatellite loci. No pattern of isolation-by-distance was detected, but a spatial principal component analysis (sPCA) revealed the presence of genetic subdivisions. Overall, the sPCA revealed two spatially distinct genetic clusters corresponding to the northern and southern parts of the study area, plus another subdivision within each of these two clusters. The north–south genetic differentiation partly matched the distribution of a snow goose colony, which could reflect a preference for settling into familiar ecological environments. Secondary clusters may result from higher-order social structures (neighbourhoods) that use landscape features to delimit their borders. The cryptic genetic subdivisions found in our population may highlight ecological processes deserving further investigations in arctic foxes at larger, regional spatial scales.

Introduction

The arctic fox Vulpes lagopus has a circumpolar distribution and inhabits all accessible Arctic tundra habitats [1]. The species is relatively common across its range, except in Fennoscandia and on islands in the Bering Sea where populations are at critically low levels [1]. At the circumpolar scale, panmixia was reported in arctic foxes, due to the combination of connectivity offered by the sea ice and high mobility of this species [2,3,4,5]. Recently, however, kin clustering linked to female philopatry and fine-scale population genetic differentiation were found in Svalbard and Alaska despite the presence of sea ice in these regions [6, 7]. Currently, most studies using microsatellites in wild arctic fox populations typically used 10 markers or less [reviewed in 8]. Genetic homogeneity of populations was found at continental or circumpolar scales, but the use of a higher genetic resolution may allow detection of fine-scale genetic structure at more regional scales [9, 10].

Building on a previous study of arctic fox extra-pair mating in the same study area [11], we developed a panel of 15 microsatellite markers combined and amplified in two multiplex and one singleplex PCR assays. We used this set of highly polymorphic microsatellite markers to assess the genetic structure of a High Arctic population in a heterogeneous landscape.

Main text

Methods

Tissue samples were collected from individuals captured on Bylot Island (73°N, 80°W), Nunavut, Canada, from 2003 to 2015, as detailed in Ref. [11]. The 600-km2 study area comprises approximately 60 km of coastline extending 5–15 km inland. The southern part of the study area hosts a Greater snow goose Chen caerulescens atlantica nesting colony during summer (Fig. 1). Total genomic DNA was extracted from ethanol-fixed ear samples stored at − 20 °C using the Qiagen DNeasy Kit. We randomly selected 15 samples from the tissue collection and 15 polymorphic primer pairs (Ref. [12,13,14,15]; Table 1) for microsatellite genotyping. These primers were used to determine optimal amplification conditions. PCRs were performed in a volume of 10 µl reaction containing 1.00 mM MgCl2, 1.00 µl DMSO, 1.00 µl BSA, 0.06 µl QBioTaq (all from Qiagen, Hilden, Germany), 0.40 µM dNTP (MP Biomedicals Europe, Illkirch, France), 0.25 µM each of forward and reverse primers (Eurofins Genomics, Ebersberg, Germany), 4.03 µl H2O and 2 ng DNA. Amplification conditions were as follows: 95 °C for 5 min, 30 cycles at 95 °C for 30 s, locus-specific annealing temperature (Ta) for 1.3 min, 72 °C for 30 s and a final extension at 60 °C for 10 min. Ta of every primer pair was determined with a Mastercycler gradient 107 thermocycler (Eppendorf) PCR machine set according to the melting temperature (Tm) values of primer pairs. The PCR products were electrophoresed on 6% denaturing polyacrylamide gels. These 15 primer pairs were selected to synthesize fluorescently-labeled primers for multiplex PCR.

Fig. 1
figure1

Map of spatial PCA global scores for arctic foxes in the south plain of Bylot Island (Nunavut, Canada) displaying a the first and b the second principal components. Each square represents the score of an individual positioned by its spatial coordinates. White squares show negative scores and black squares show positive scores. Larger squares reflect greater absolute values. Squares of different colors are strongly differentiated, while squares of the same color but of different sizes are weakly differentiated. The study area is delimited by a black line. Blue lines indicate rivers. The extent of a snow goose nesting colony is shown in red. Base map source from Google Maps: Imagery ©2017 Google, TerraMetrics

Table 1 Summary information for the 15 microsatellite markers used in this study of Vulpes lagopus

A Qiagen multiplex PCR Kit was used and reaction mixtures contained 1 µl DNA, 6.25 µl Master Mix (Qiagen, Hilden, Germany) with 1.25 µl primer mix (100 pM/µl) and 3.5 µl RNAse-free water to a final volume of 12 µl. The reaction comprised an activation step at 95 °C for 5 min, followed by 35 cycles of initial denaturation at 95 °C for 30 s, Ta for 30 s (Table 1) and 72 °C for 30 s, ending with a final extension step at 60 °C for 30 min, and followed by a holding step at 20 °C. Fluorescently-labeled PCR products were run on an ABI PRISM 3130 DNA Sequencer (Applied Biosystems) with the GeneScan-500 (LIZ) internal size standard and analyzed with the GeneMapper software (Applied Biosystems).

Fifteen loci (Multiplex1: 9 loci; Multiplex2: 5 loci; Singleplex: 1 locus; Table 1) were successfully amplified and genotyped in 210 specimens. We used R package adegenet_2.0.1 [16] to estimate the number of alleles (NA), gene diversity He (expected heterozygosity), and observed heterozygosity (Ho) at each locus (Table 1). We used GENEPOP 4.1.3 [17] with a Markov chain method (10,000 dememorization steps, 100 batches and 5000 iterations per batch) to calculate inbreeding coefficients (FIS) and test for Hardy–Weinberg equilibrium (HWE) and linkage disequilibrium. Significance levels were adjusted for multiple tests using the sequential Bonferroni technique [18]. We used ML-NULLFREQ to estimate null allele frequencies based on maximum likelihood methods [19].

We investigated the fine-scale genetic structure using a global spatial autocorrelation technique (Mantel correlograms) [20, 21] implemented in the software GenAlEx 6.5 [22]. We used data of genotyped individuals that were resident in the study area (seven individuals dispersed, therefore n = 203). Spatial coordinates corresponded to occupied dens or capture sites. We used even distance classes of 4 km since this distance represents the mean radius of an arctic fox’s home range [23]. Statistical testing was based on 999 random permutations of individual genotypes. The complete dataset was first used, and then divided by sex (n = 99 males and 104 females) to test for sex-biased kin clustering. In addition, we used a spatially explicit multivariate method, the spatial principal component analysis (sPCA) implemented in adegenet_2.0.1. The sPCA takes into account both the genetic diversity based on allele frequencies (variance) and the spatial autocorrelation between individuals (Moran’s I). Global structures, such as distinct patches, genetic clines or intermediate states, appear when immediate neighbors are more genetically similar than expected under a random distribution. Local structures detect strong genetic differences among immediate neighbors, which could arise when individuals with a similar genetic pool avoid mating with each other. This approach does not require genetic equilibrium conditions and is useful for revealing cryptic spatial patterns [24], as demonstrated in other mammalian carnivores [20, 21, 25, 26]. We used a distance-based connection network, with a threshold distance between any two neighbors chosen as the minimum distance so that no individual was excluded from the graph [24]. As individuals could be associated to the same den, we added 100 m (in a random direction) to individual coordinates. We tested for significant global and local structure using a Monte-Carlo randomization test (999 permutations). We used the complete dataset first (threshold distance = 4359 m), and then the sex-specific datasets (threshold distance = 4766 m for males and 4424 m for females). For comparison, we also performed a Bayesian cluster analysis using the software STRUCTURE 2.3.4 [27].

Results and discussion

The number of alleles per locus varied from five to 17. Expected heterozygosity ranged from 0.667 to 0.890, while observed heterozygosity ranged from 0.576 to 0.886. Null alleles occurred at low frequency (0.000–0.020) for 12 loci and at medium frequency for three loci (0.038–0.090). Significant departures from Hardy–Weinberg equilibrium after Bonferroni correction were detected at four loci (Table 1). Eight pairs of loci (7.6%) deviated significantly from linkage equilibrium after Bonferroni correction, a number only slightly higher than what may be expected by chance alone (5.25). As these deviations were very likely the result of the significant spatial genetic structure detected in the subsequent analyses, we included all loci in the analyses.

No significant spatial genetic structure was detected by the Mantel correlogram using the complete dataset (Mantel r = − 0.023 to 0.003; all p > 0.066) or sex-specific datasets (males: Mantel r = − 0.032 to 0.016, all p > 0.097; females: Mantel r = − 0.026 to 0.005; all p > 0.112). These results contrast with those of Ehrich et al. [6], who reported an isolation-by-distance pattern and female kin clustering in Svalbard. The observed trend was however very weak in their study (R2 = 0.0007, p = 0.027) and the median distance reported between first-order relatives (35.7 km) was relatively high [6].

The Monte-Carlo tests from the sPCA using the complete dataset revealed the existence of at least one global pattern (observed value = 0.009, p = 0.003) but no local structure (observed value = 0.010, p = 0.445). The first two eigenvalues were retained after examination of the eigenvalues barplot and the sPCA screeplot (Additional file 1: Figure S1). The scores from the first principal component mainly showed a latitudinal differentiation in the study area, with the southern cluster matching closely the extent of the goose colony (Fig. 1a). The scores from the second principal component mainly separated the individuals in the northernmost and southernmost parts of the study area from the rest, with the two clusters in the northern part (north of the goose colony) showing a sharper boundary (larger squares) than in the south (Fig. 1b). The southern part of the study area (goose colony and south of the colony) presented a north–south cline (genotypes located in the middle of the distribution had less extreme scores, as depicted by the smaller squares), indicating more mixing between the individuals in southern part of the goose colony.

Such cryptic population genetic structure in areas where no distinct physical barriers to movement occur may reveal interesting ecological processes. For instance, the principal southern cluster which approximately matched the goose nesting colony may reflect habitat imprinting (also termed, natal habitat-biased dispersal), which is a preference for settling into familiar ecological environments [28, 29]. This behavior has been suggested for arctic foxes to explain why the “lemming ecotype”, which feeds preferentially on small rodents, seems less likely to settle in areas inhabited by the “coastal ecotype”, which relies on resources from sea bird colonies and the marine environment [9, 30]. Habitat imprinting based on prey types or habitat features leading to cryptic genetic subdivisions has been reported in other canids, such as coyotes Canis latrans [31, 32] and gray wolves Canis lupus [21]. The smaller subdivisions revealed by the second principal component axis may be more difficult to explain. Their boundaries coincide approximately with rivers bordered by steep hillsides (Fig. 1b), but it is unlikely that these imped fox movements as foxes are regularly observed swimming across rivers (personal observations). However, the most distinct clusters may represent higher-order social structures (neighbourhoods occupied by extended family members) that may have borders delimited by landscape features such as rivers. Such an interaction between landscape features and social cohesion may become apparent at fine scales, as seen in the central California coyote population [32]. While remaining territorial, breeding pairs may adopt a good-neighbor strategy by tolerating relatives as neighbours, as observed in Norway where close relatives sometimes engage in communal breeding [33]. The sPCAs with either males or females did not reveal any pattern (data not shown), indicating no sex-biased genetic structures. The Bayesian cluster analysis inferred K = 2 as the most likely number of clusters (Additional file 2: Figure S2a). The percentage of individuals with membership coefficients q > 0.7 was relatively low (61%), highlighting a high level of admixture. Overall, individuals assigned to Cluster I and with a mixed membership were distributed all over the study area, while 73% of those assigned to Cluster II were found in the goose colony area (Additional file 2: Figure S2b, c).

The influence of social relationships on natal dispersal and spacing patterns, however, requires further investigation. From 2003 to 2015, a total of 371 pups were tagged in the population [annual pup production tagged (mean ± SD): 47 ± 27%]. Fourteen of these were resighted as adults and only six bred in the study area. Although not all pups in the study area were tagged, this suggests a lack of natal philopatry (or high mortality).

Conclusions

This study confirms the absence of an isolation-by-distance pattern but reveals cryptic genetic subdivisions in our population. As reported in other studies [21, 31], using linear distances to estimate canid dispersal patterns may be misleading, especially in heterogeneous landscapes, and spatially explicit multivariate methods such as the sPCA may be more appropriate. Our multiplex assay consisting of 15 polymorphic loci will enhance analytical power to conduct fine-scale genetic analyses for this highly studied and locally endangered species which is the focus of > 30 monitoring projects throughout the world [34]. A comprehensive understanding of the spatial distribution of genetic diversity will inform conservation efforts.

Limitations

One limitation of our study is that the study area may have been too small to detect a clear pattern of isolation-by-distance for a mobile species. Considering the relatively large distances between close relatives reported in Svalbard [6], foxes in sea ice regions may not display strict philopatry (staying in or moving to a site adjacent to the natal home range) as it is observed in Scandinavia [35]. Although logistically difficult to perform, increasing the sampling area may allow us to capture philopatry at a larger scale.

References

  1. 1.

    Angerbjörn A, Hersteinsson P, Tannerfeldt M. Arctic fox (Alopex lagopus). In: MacDonald DW, Sillero-Zubiri C, editors. Canids: foxes, wolves, jackals and dogs—2004 status survey and conservation action plan. Gland: IUCN; 2004. p. 117–23.

  2. 2.

    Norén K, Carmichael LE, Dalén L, Hersteinsson P, Samelius G, Fuglei E, Kapel CMO, Menyushina I, Strobeck C, Angerbjörn A. Arctic fox Vulpes lagopus population structure: circumpolar patterns and processes. Oikos. 2011;120:873–85.

  3. 3.

    Carmichael LE. Genetics of North American arctic fox populations. In: Final wildlife report, vol. 16. Iqaluit: Government of Nunavut, Department of Environment; 2006. p. 1–37.

  4. 4.

    Dalen 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.

  5. 5.

    Geffen E, Waidyaratne S, Dalen L, Angerbjörn A, Vila C, Hersteinsson P, Fuglei E, White PA, Goltsman M, Kapel CMO, Wayne RK. Sea ice occurrence predicts genetic isolation in the Arctic fox. Mol Ecol. 2007;16:4241–55.

  6. 6.

    Ehrich D, Carmichael L, Fuglei E. Age-dependent genetic structure of arctic foxes in Svalbard. Polar Biol. 2012;35:53–62.

  7. 7.

    Goldsmith EW, Renshaw B, Clement CJ, Himschoot EA, Hundertmark KJ, Hueffer K. Population structure of two rabies hosts relative to the known distribution of rabies virus variants in Alaska. Mol Ecol. 2016;25:675–88.

  8. 8.

    Norén K, Dalén L, Flagstad Ø, Berteaux D, Wallén J, Angerbjörn A. Evolution, ecology and conservation—revisiting three decades of arctic fox population genetic research. Polar Res. 2017;36:4.

  9. 9.

    Norén K, Carmichael L, Fuglei E, Eide NE, Hersteinsson P, Angerbjörn A. Pulses of movement across the sea ice: population connectivity and temporal genetic structure in the arctic fox. Oecologia. 2011;166:973–84.

  10. 10.

    Meinke PG, Kapel CMO, Arctander P. Genetic differentiation of populations of Greenlandic Arctic fox. Polar Res. 2001;20:75–83.

  11. 11.

    Cameron C, Berteaux D, Dufresne F. Spatial variation in food availability predicts extrapair paternity in the arctic fox. Behav Ecol. 2011;22:1364–73.

  12. 12.

    Fredholm M, Winterø AK. Variation of short tandem repeats within and between species belonging to the Canidae family. Mamm Genome. 1995;6:11–8.

  13. 13.

    Ostrander EA, Sprague GF Jr, Rine J. Identification and characterization of dinucleotide repeat (CA)n markers for genetic mapping in dog. Genomics. 1993;16:207–13.

  14. 14.

    Ostrander EA, Mapa FA, Yee M, Rine J. One hundred and one simple sequence repeat-based markers for the canine genome. Mamm Genome. 1995;6:192–5.

  15. 15.

    Mellersh CS, Langston AA, Acland GM, Fleming MA, Ray K, Wiegand NA, Francisco LV, Gibbs M, Aguirre GD, Ostrander EA. A linkage map of the canine genome. Genomics. 1997;46:326–36.

  16. 16.

    Jombart T. adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24:1403–5.

  17. 17.

    Rousset F. Genepop’007: a complete re-implementation of the genepop software for Windows and Linux. Mol Ecol Resour. 2008;8:103–6.

  18. 18.

    Rice WR. Analyzing tables of statistical tests. Evolution. 1989;43:223–5.

  19. 19.

    Kalinowski ST, Taper ML. Maximum likelihood estimation of the frequency of null alleles at microsatellite loci. Conserv Genet. 2006;7:991–5.

  20. 20.

    Reding DM, Bronikowski AM, Johnson WE, Clark WR. Pleistocene and ecological effects on continental-scale genetic differentiation in the bobcat (Lynx rufus). Mol Ecol. 2012;21:3078–93.

  21. 21.

    Stronen AV, Navid EL, Quinn MS, Paquet PC, Bryan HM, Darimont CT. Population genetic structure of gray wolves (Canis lupus) in a marine archipelago suggests island-mainland differentiation consistent with dietary niche. BMC Ecol. 2014;14:11.

  22. 22.

    Peakall R, Smouse PE. GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research—an update. Bioinformatics. 2012;28:2537–9.

  23. 23.

    Lai S, Bêty J, Berteaux D. Movement tactics of a mobile predator in a meta-ecosystem with fluctuating resources: the arctic fox in the High Arctic. Oikos. 2017;126:937–47.

  24. 24.

    Jombart T, Devillard S, Dufour AB, Pontier D. Revealing cryptic spatial patterns in genetic variability by a new multivariate method. Heredity. 2008;101:92–103.

  25. 25.

    Basto MP, Santos-Reis M, Simões L, Grilo C, Cardoso L, Cortes H, Bruford MW, Fernandes C. Assessing genetic structure in common but ecologically distinct carnivores: the stone marten and red fox. PLoS ONE. 2016;11:e0145165.

  26. 26.

    Warren MJ, Wallin DO, Beausoleil RA, Warheit KI. Forest cover mediates genetic connectivity of northwestern cougars. Conserv Genet. 2016;17:1011–24.

  27. 27.

    Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.

  28. 28.

    Wecker SC. The role of early experience in habitat selection by the prairie deer mouse, Peromyscus maniculatus bairdi. Ecol Monogr. 1963;33:307–25.

  29. 29.

    Vogl W, Taborsky M, Taborsky B, Teuschl Y, Honza M. Cuckoo females preferentially use specific habitats when searching for host nests. Anim Behav. 2002;64:843–50.

  30. 30.

    Pagh S, Hersteinsson P. Difference in diet and age structure of blue and white Arctic foxes (Vulpes lagopus) in the Disko Bay area, West Greenland. Polar Res. 2008;27:44–51.

  31. 31.

    Sacks BN, Brown SK, Ernest HB. Population structure of California coyotes corresponds to habitat specific breaks and illuminates species history. Mol Ecol. 2004;13:1265–75.

  32. 32.

    Sacks BN, Mitchell BR, Williams CL, Ernest HB. Coyote movements and social structure along a cryptic population genetic subdivision. Mol Ecol. 2005;14:1241–9.

  33. 33.

    Elmhagen B, Hersteinsson P, Norén K, Unnsteinsdottir E, Angerbjörn A. From breeding pairs to fox towns: the social organisation of arctic fox populations with stable and fluctuating availability of food. Polar Biol. 2014;37:111–22.

  34. 34.

    Berteaux D, Thierry A-M, Alisauskas R, Angerbjörn A, Buchel E, Doronina L, Ehrich D, Eide NE, Erlandsson R, Flagstad Ø, et al. Harmonizing circumpolar monitoring of Arctic fox: benefits, opportunities, challenges and recommendations. Polar Res. 2017;36:2.

  35. 35.

    Tannerfeldt M, Angerbjörn A. Life history strategies in a fluctuating environment: establishment and reproductive success in the arctic fox. Ecography. 1996;19:209–20.

Download references

Authors’ contributions

SL, DB and AL participated in the design of the study; SL and DB planned the field work; SL and AL carried out data collection; AL, JL and AQ identified polymorphic loci and genotyped the samples; SL and AL performed the analyses and wrote the manuscript. All authors read and approved the final manuscript.

Acknowledgements

We thank our many field assistants and students who helped collect the data and Justin Roy (UQAR) for advice on the analyses. We also thank the ‘Service de Systématique Moléculaire’ (UMS2700 MNHN/CNRS) for granting access to its technical platform.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

DNA sequences of the microsatellite markers are included in the article. The raw data generated and analysed in this study are available from the corresponding author upon request.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Fox capture and handling procedures were approved by the Animal Care Committee of Université du Québec à Rimouski [permit # CPA32-08-62(R2)].

Funding

This study was supported by Canada Foundation for Innovation, Canada Research Chairs Program, Fonds de recherche du Québec-Nature et technologies (FRQNT), International Polar Year program, Kenneth M Molson Foundation, Mittimatalik Hunters and Trappers Organization, Natural Sciences and Engineering Research Council of Canada (NSERC), Network of Centers of Excellence of Canada ArcticNet, Northern Scientific Training Program (Polar Knowledge Canada), Nunavut Wildlife Management Board, Parks Canada Agency, Polar Continental Shelf Program (Natural Resources Canada), Université du Québec à Rimouski (UQAR), and the ‘UMR ISYEB MNHN grant 2014’.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information

Correspondence to Sandra Lai.

Additional files

13104_2017_3002_MOESM1_ESM.pdf

Additional file 1: Figure S1. Barplot of spatial principal component analysis (sPCA) eigenvalues and screeplot displaying each eigenvalue according to its variance and spatial autocorrelation (Moran’s I) components. Graphical results of the sPCA analysis of arctic foxes (n = 203) from Bylot Island, Nunavut, Canada.

13104_2017_3002_MOESM2_ESM.pdf

Additional file 2: Figure S2. Bayesian genetic structure analysis conducted with the software STRUCTURE. Graphical results of the Bayesian genetic structure analysis of arctic foxes (n = 203) from Bylot Island, Nunavut, Canada conducted with the software STRUCTURE.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • Vulpes lagopus
  • Microsatellite multiplex PCR
  • Population genetics
  • Fine-scale genetic structure
  • Bylot Island