Museomics Dissects the Genetic Basis for Adaptive Seasonal Coloration in the Least Weasel.

Dissecting the link between genetic variation and adaptive phenotypes provides outstanding opportunities to understand fundamental evolutionary processes. Here, we use a museomics approach to investigate the genetic basis and evolution of winter coat coloration morphs in least weasels (Mustela nivalis), a repeated adaptation for camouflage in mammals with seasonal pelage color moults across regions with varying winter snow. Whole-genome sequence data were obtained from biological collections and mapped onto a newly assembled reference genome for the species. Sampling represented two replicate transition zones between nivalis and vulgaris coloration morphs in Europe, which typically develop white or brown winter coats, respectively. Population analyses showed that the morph distribution across transition zones is not a by-product of historical structure. Association scans linked a 200-kb genomic region to coloration morph, which was validated by genotyping museum specimens from intermorph experimental crosses. Genotyping the wild populations narrowed down the association to pigmentation gene MC1R and pinpointed a candidate amino acid change cosegregating with coloration morph. This polymorphism replaces an ancestral leucine residue by lysine at the start of the first extracellular loop of the protein in the vulgaris morph. A selective sweep signature overlapped the association region in vulgaris, suggesting that past adaptation favored winter-brown morphs and can anchor future adaptive responses to decreasing winter snow. Using biological collections as valuable resources to study natural adaptations, our study showed a new evolutionary route generating winter color variation in mammals and that seasonal camouflage can be modulated by changes at single key genes.

Introduction

Dissecting the relationship between adaptive phenotypes and the underlying genetic variation is crucial to understand how adaptations emerge and are maintained across populations and the intricate relationship between the distribution of trait variation and environmental selective pressures (Savolainen et al. 2013). Seasonal color change is a key phenological adaptation of over 20 mammal and bird species (Mills et al. 2018). Seasonal moults between summer-brown and winter-white coats, mainly regulated by photoperiod, allow camouflage year-round in environments periodically covered with snow (Mills et al. 2013; Zimova et al. 2018). Within these species, genetically determined polymorphism in winter color, that is, winter-white or winter-brown morphs, often varies clinically across the natural distribution range of the species and correlates with snow cover variables (Mills et al. 2018; Zimova et al. 2018). The fitness cost of coat-background color mismatches (Imperio et al. 2013; Sultaire et al. 2016; Zimova et al. 2016; Atmeh et al. 2018; Wilson et al. 2018; Melin et al. 2020) can be particularly harmful to winter-white populations facing decreased snow cover duration driven by global warming (Mills et al. 2013). Identifying the molecular mechanisms underlying polymorphism in seasonal coat color is therefore crucial to understand not only the evolutionary pathways for the repeated emergence of adaptive trait variation (Zimova et al. 2018; Giska et al. 2019) but also the potential for rapid adaptation in response to climate change (Mills et al. 2018). Yet, our understanding of the evolution of this adaptive polymorphism across species and genera remains limited (Våge et al. 2005; Jones et al. 2018; Giska et al. 2019; Jones et al. 2020a). In this work, we dissect the genetic basis of adaptive color polymorphism in the least weasel (Mustela nivalis).

Least weasels are small carnivores with a wide Holarctic distribution (Abramov and Baryshnikov 2000). In Europe, two subspecies are consistently described, M. n. nivalis and M. n. vulgaris (Abramov and Baryshnikov 2000; King and Powell 2007), but this taxonomic classification is not coherent with population structure inferred from mitochondrial DNA variation, in particular at transition zones between morphs (Lebarbenchon et al. 2010; McDevitt et al. 2012). The subspecies are rather distinguished by morphological traits: size and coloration type ( fig. 1 and supplementary fig. S1, Supplementary Material online). The nivalis morph is characterized by a typical winter-white coat, a straight demarcation line separating the dorsal and ventral pelages, and smaller size. The vulgaris morph does not develop winter-white coats, presents a ragged dorsoventral line, and is larger (Abramov and Baryshnikov 2000; King and Powell 2007). Since the terms nivalis and vulgaris are commonly used for coloration morphs, also designated as type I and type II, respectively (Frank 1985), we will do so hereafter. Experimental crosses between the nivalis and vulgaris morphs did not separate the segregation of the coloration traits (dorsoventral line and color moult), which were suggested to be determined by the same locus or by tightly linked loci (Frank 1985). This also showed that coloration morphs have Mendelian inheritance, with vulgaris being dominant and nivalis recessive (Frank 1985). In Europe, the nivalis morph is found at higher latitudes and altitudes where winter snow cover is more persistent, whereas the vulgaris morph inhabits southern and coastal areas. This clinal morph distribution correlates with snow cover-related variables (Mills et al. 2018), suggesting that the distribution of the coloration types is maintained by environmental selection for camouflage. In the transition zones between these distributions, the coloration morphs coexist (Mills et al. 2018; fig. 1).

Fig. 1.

Coat color polymorphism and population structure in the least weasel. (A) Mustela nivalis individuals in (from left to right) typical summer nivalis coloration, winter nivalis coloration, and year-round vulgaris coloration (Photos: Karol Zub). (B) Clinal distribution of coloration morphs shown as the probability of winter-white coats across the European distribution range of the species (adapted from Mills et al. 2018, with data available from Mills et al. 2019). Striped shading depicts areas outside the species’ natural distribution defined in the IUCN database (http://www.iucnredlist.org, last accessed April 11, 2020). (C) Zoom-in of the sampling area, delimited by the dashed line in (B), indicating the localities of analyzed specimens in Sweden and Poland (overlapping sample localities are indicated by a single symbol). Circles denote specimens analyzed with whole-genome sequencing and genotyping data, whereas diamonds denote specimens used only for genotyping. Symbol colors indicate sampled coloration morphs (brown—vulgaris, white—nivalis, brown/white—both morphs). (D) Admixture proportions inferred with K = 2 genetic clusters based on 71,867 unlinked SNPs. Coloration morphs are shown as white (nivalis) or brown (vulgaris) circles. (E) Population/morph tree derived from allele frequencies estimated for 118,642 unlinked genome-wide SNPs. The ferret (M. putorius furo) was used as outgroup. All nodes have bootstrap support of 1.

Our modern ability to apply high-throughput DNA-sequencing technologies to historical samples provides unprecedented precision to infer biological processes from museum data (Burrell et al. 2015) and to sample adaptive phenotypic variation. Natural history collections are rich repositories of biological variation at inter- and intra-specific level, which provide fundamental opportunities to study evolutionary processes at multiple taxonomic scales (Holmes et al. 2016; Meineke et al. 2019), allowing, for example, molecular-based biodiversity assessments (Mikheyev et al. 2017), phylogenetic reconstructions (Wood et al. 2018), or temporal analyses of genetic variation (Bi et al. 2013, 2019; van der Valk et al. 2019). Here, we expand the potential of museomics and use whole-genome scans in specimens from natural history collections to identify the genetic basis and assess the role of natural selection in the evolution of the nivalis and vulgaris coloration morphs, in wild populations of the least weasel. Our results provide important insights into the repeated evolution of a remarkable adaptive polymorphism and the selective mechanisms maintaining geographical clines of adaptive variants. Further, we underscore the importance of biological collections as crucial resources to study adaptive phenotypes.

Results and Discussion

Population Structure Is Determined by Geography Rather than Coloration Morph

Whole genomes from 79 least weasel specimens sampled across two coloration morph transition zones in Sweden (N = 20 nivalis, N = 20 vulgaris) and Poland (N = 20 nivalis, N = 19 vulgaris) (Atmeh et al. 2018; Mills et al. 2018; see Materials and Methods for a detailed characterization of the phenotypic classes; fig. 1C and supplementary table S1, Supplementary Material online) were sequenced at low individual coverage (1.05× mean coverage in Swedish specimens, 1.38× in Polish specimens; supplementary table S2, Supplementary Material online) from individually barcoded DNA libraries. The sampled specimens were retrieved from the collections of the Swedish Museum of Natural History (NRM) and the Mammal Research Institute of the Polish Academy of Sciences (MRI) and were originally collected between 1959 and 2010 (supplementary table S1, Supplementary Material online). The nivalis/vulgaris winter coloration phenotype was recorded from the prepared skins of the Swedish specimens, whereas for the Polish samples, the nivalis/vulgaris morph was recovered from the original data record and the geographical locality of the sample, considering the known distribution of the morphs (Atmeh et al. 2018). Illumina reads were mapped onto a de novo assembly of the least weasel genome, which we reconstructed from a male road-killed nivalis specimen opportunistically sampled in Poland, using Chromium 10× linked-read data (length 2.50 Gb; N50 = 24.90 Mb; 33,286 scaffolds; 95.3% completeness based on 4,104 core mammalian genes; see supplementary table S3, Supplementary Material online). Principal component analysis (PCA), individual clustering inferences, and permutational analyses of variance between morphs based on 71,867 unlinked genome-wide single-nucleotide polymorphisms (SNPs) showed that geography is the major determinant of population structure, separating primarily the Swedish from the Polish specimens, regardless of coloration morph (PERMANOVA, P > 0.05; fig. 1 and supplementary fig. S2A, Supplementary Material online). Similarly, clustering followed geography rather than coloration morph in a neighbor-joining tree inferred from individual pairwise genetic distances (supplementary fig. S2B, Supplementary Material online), in a population tree inferred from genome-wide allele counts (118,642 sites; fig. 1), and in a complete mitochondrial DNA phylogeny (31 specimens from Sweden, 33 from Poland; supplementary fig. S3 and supplementary table S2, Supplementary Material online). These results suggest that the overall population structure does not reflect coloration morph distributions and thus appears unrelated to the subspecific classification. At the local scale, shallow structure was found in both Swedish and Polish populations (93,374 and 87,803 SNPs, respectively; supplementary fig. S4, Supplementary Material online), with no principal component significantly explaining the variance of the genetic data (randomization test, P > 0.05). In both localities, mean genome-wide FST between the nivalis and vulgaris morphs was low (∼0.068 in Sweden based on 7,330,638 SNPs and ∼0.059 in Poland based on 10,813,464 SNPs), confirming the shallow structure. However, permutational analyses of variance suggested that the centroids of the genome-wide variation differ between morphs in Sweden and Poland (PERMANOVA; P < 0.001 and P < 0.01, respectively), which may suggest some degree of historical differentiation between morphs at the local scale. Regardless of the local specificities of environmental-related transitions in the distribution of coloration morphs (Mills et al. 2018), our analyses show that the Swedish and Polish populations of least weasels are two replicates of the nivalis–vulgaris phenotypic transition, providing a particularly suitable model to map the genetic basis of the coloration trait.

MC1R Variation Is Associated with Coloration Morph

Taking advantage of the two replicated transition zones sampled, we performed a genome-wide scan of the Cochran–Mantel–Haenszel (CMH) association test (based on 7,821,284 SNPs inferred from pooled data; McDonald 2009) and identified consistent differences in allele frequencies between coloration morphs in Sweden and Poland clustered in scaffold 337 (fig. 2). Similar results were obtained from a Fisher exact test (FET) of allele frequency differences (13,268,106 SNPs; supplementary fig. S5A, Supplementary Material online) and a case–control association test (17,685,066 SNPs estimated from genotype likelihoods; fig. 2), and the association remained when controlling for global population structure in an association score test (4,638,022 SNPs; supplementary fig. S5B, Supplementary Material online). This outlier genomic region was also recovered when analyzing the Swedish and Polish transitions separately (FET and case-control association scans; supplementary fig. S5C–F, Supplementary Material online). Window-based FST scans between morphs showed a longer region of elevated differentiation in scaffold 337 in Sweden (fig. 2) than in Poland (fig. 2 and supplementary fig. S5G, Supplementary Material online), suggesting a longer stretch of linkage disequilibrium across the association region in the Swedish population. Taken together, these results point to a genomic region of ∼200 kb of elevated differentiation in scaffold 337 (coordinates 125,000–325,000 bp) as possibly underlying the nivalis–vulgaris coloration differences in least weasels. This region is homologous to a region in scaffold GL896939.1 of the ferret genome (between positions 134,432–337,301 bp; Peng et al. 2014), which includes seven genes: DEF8, TUBB3, MC1R, TCF25, SPIRE2, FANCA, and ZNF276 (fig. 2). The melanocortin-1 receptor (MC1R) is a well-known pigmentation gene (García-Borrón et al. 2005), involved in numerous cases of color polymorphism in wild and domesticated organisms (Hubbard et al. 2010). The spire-type actin nucleation factor 2 (SPIRE2) gene has been shown to be involved in the regulation of myosin-Va-dependent melanosome dispersion in melanocytes (Alzahofi et al. 2020), which could indicate an indirect function in pigment distribution. None of the other genes has a known function in coloration.

Fig. 2.

Genetic mapping of the association with coloration morph in the least weasel. (A) CMH test of significance of allele frequency differences between coloration morphs across the Swedish and Polish populations of the least weasel. The red line indicates the Bonferroni-corrected threshold of P = 0.05. (B) Case–control association test of allele frequency differences between nivalis and vulgaris for the complete data set (Swedish and Polish populations). The red line indicates the Bonferroni-corrected threshold of P = 0.05. (C) Whole-genome scan of FST between nivalis and vulgaris coloration morphs in Sweden. Values are averaged across 25 kb nonoverlapping windows. The red line represents the 99.9th percentile of FST values across windows. In panels (A–C), scaffolds were ordered from left to right by decreasing length. (D) Zoom-in of FST values between coloration morphs for the Swedish (circles, averaged in 25 kb windows with 5 kb steps) and Polish (diamonds, averaged in 5 kb nonoverlapping windows) populations in the first 500 kb of scaffold 337. Orange circles and blue diamonds highlight windows in the 99th FST percentile in the Swedish and Polish populations, respectively. Genes annotated along the region are shown as rectangles and those within the candidate region are identified. (E) Genotypes at 40 loci spanning the candidate region including MC1R (scaffold 337) and two additional outliers of differentiation (scaffolds 333 and 338). Each row depicts a specimen and each column a genotyped locus. The left column indicates the population of origin and coloration morphs (white—nivalis or brown—vulgaris). For scaffold 337, the distribution of SNPs across genes is shown. The dashed black box highlights the candidate SNPs identified in the MC1R exon. For each population, genotypes are colored as light gray—homozygous for the predominant nivalis variant; black—homozygous for the alternative variant; dark gray—heterozygous; white—missing data.

To validate this association, we used an independent sample set composed of 35 individuals from the interbreeding experiment between the nivalis and vulgaris morphs performed by Frank (1985), deposited at the collection of the Zoological Research Museum Alexander Koenig (Bonn, Germany; supplementary table S4, Supplementary Material online). Frank (1985) crossed a nivalis female imported from Sweden with a vulgaris male from Germany, and backcrossed F1 and BC1 (first backcross) male offspring with the nivalis female, recording the individual coloration morphs. We genotyped this sample set for 22 SNPs showing elevated allele frequency differences (mainly >0.5) between morphs in the Swedish transition zone, covering the longest association tract found in scaffold 337 (of which 19 were also polymorphic in the Polish population). Eighteen additional SNPs located in other windows of high differentiation between morphs found in our genome scans for the Swedish transition (top 0.01% of FST values; scaffolds 333 and 338; fig. 2) were also genotyped. Frank (1985) inferred that the coloration morph has Mendelian segregation, being nivalis recessive and vulgaris dominant. Therefore, all nivalis specimens are expected to be homozygous at causal or linked loci, whereas vulgaris may either be heterozygous or homozygous for the alternative variant. Eighteen out of the 22 SNPs from scaffold 337 followed the predicted genotypes across all analyzed specimens from Frank (1985), as expected from a large linkage block resulting from early generation crosses, confirming the association with coloration morph (vulgaris dominant model; P < 0.05 corrected for multiple tests; fig. 2 and supplementary tables S5 and S6, Supplementary Material online). In contrast, the additional SNPs located in scaffolds 333 and 338 did not (P > 0.05 correct for multiple tests; fig. 2 and supplementary tables S5 and S6, Supplementary Material online), suggesting that the high differentiation found in the Swedish population in these other genomic regions was likely a by-product of local population structure unrelated with coloration differences.

To narrow down the association, we took advantage of the expected reduction in linkage disequilibrium across the region of association in the wild and genotyped the same SNP set in extended population samples (N = 80 in Sweden, NRM collection; N = 58 in Poland, MRI collection; supplementary table S1, Supplementary Material online). All 22 SNPs located across the association region in scaffold 337 were found significantly associated with coloration morph in Sweden and ten in Poland (P < 0.05 corrected for multiple tests; fig. 2 and supplementary table S6, Supplementary Material online). Among the 25 Polish vulgaris specimens genotyped, we found that 12 share a haplotype block similar to Swedish vulgaris, whereas in 13, the association block is shortened (fig. 2 and supplementary table S5, Supplementary Material online), which explains the shorter association region identified from the genome scans (fig. 2). We inspected the genotyped SNPs along the association region and found that only two consecutive SNPs follow the inheritance pattern predicted from Frank (1985) (fig. 2 and supplementary table S5, Supplementary Material online), not only in specimens for which we could score the coloration morph but also those for which we assumed an a priori nivalis/vulgaris identification recorded in the collection database. These two variants lay on the MC1R reading frame (supplementary table S7, Supplementary Material online), suggesting that this gene plays a role in the coloration differences studied here.

The MC1R gene encodes the melanocortin-1 receptor (MC1R), a G-protein coupled transmembrane receptor of melanocytes that plays a key role in melanogenesis (García-Borrón et al. 2005). The α-melanocyte-stimulating hormone (α-MSH) is an agonist of MC1R (Rouzaud et al. 2003) and triggers the production of dark eumelanin (brown to black pigment) during melanogenesis, whereas the binding of the agouti signaling protein (ASIP) (Lu et al. 1994) shifts melanogenesis to the production of lighter pheomelanin (yellow to red pigment) or inhibits pigment production (Le Pape et al. 2009). Whole-genome analyses in hares identified ASIP as the gene causing winter coat color polymorphism in two species (white-brown and white-gray variation; Jones et al. 2018; Giska et al. 2019). Additionally, a candidate gene approach in arctic foxes proposed two amino acid changes in MC1R as underlying the blue morph, which, in contrast with the least weasel, shows clear differences in both winter and summer coats compared with the wild type (Våge et al. 2005). The ASIP–MC1R system therefore appears implicated in determining seasonal coloration differences across mammals. In hares, these differences have been shown to be determined by cis-regulatory changes influencing the expression of ASIP during the autumn moult (causing higher expression in the winter-white morphs; Jones et al. 2018; Giska et al. 2019). We further used our genomic data set to inspect the genetic variation flanking the mutations showing perfect association with coloration morph and found no other variants with increased allele frequency differences between morphs in both transition zones, including in the regulatory region of MC1R (supplementary fig. S6, Supplementary Material online). Although we cannot completely exclude the existence of undetected variation, our evidence suggests that the coloration polymorphism in M. nivalis, contrary to what was found in hares (Jones et al. 2018; Giska et al. 2019), is likely not driven by cis-regulatory changes. Understanding whether MC1R regulation may be involved in winter color polymorphism must, however, await gene expression analyses (e.g., Jones et al. 2018), which would require targeted de novo sampling of fresh tissues.

We then focused on the mutations showing cosegregation with coloration type. These two consecutive nucleotide changes in the reading frame of MC1R cause an amino acid change between leucine and lysine at residue 101 (nivalis L101K vulgaris; positions 165,924 and 165,925 in scaffold 337; supplementary fig. S7, Supplementary Material online). Both nucleotide substitutions are required for this amino acid change, and the residue was found to be conserved at a deep evolutionary scale (based on 100 species, including the major mammal groups; supplementary fig. S7A and B and supplementary table S8, Supplementary Material online), with leucine (nivalis) being the ancestral state. Only one additional amino acid with similar biochemical properties (Betts and Russell 2003)—valine—is found at this position in dogs and foxes (supplementary fig. S7A, Supplementary Material online). Tests of functional impact, based on sequence conservation and amino acid properties, suggested that the leucine-lysine change at position 101 impacts protein function (SIFT score = 0.00). The L101K substitution in the least weasel is located at the transition between the second transmembrane domain (TMD2) and the first extracellular loop of the receptor (supplementary fig. S7C, Supplementary Material online), considering the 2D structure of MC1R inferred from the human protein (Ringholm et al. 2004). Mutations occurring within or near TMD2 and TMD3 have been linked with the expression of darker (eumelanized) phenotypes in mammals (Robbins et al. 1993; Klungland et al. 1995; Våge et al. 1997; Kijas et al. 1998; Våge et al. 1999; Newton et al. 2000) and constitutive decreased efficiency of ligand binding (Lu et al. 1998). Although a constitutive activation of MC1R appears unlikely, decreased binding efficiency of the agouti signaling protein to vulgaris MC1R could elevate the threshold for the antagonist-induced response, impeding the response to the ASIP pulse of expression during the autumn moult (Ferreira et al. 2017; Jones et al. 2018; Ferreira et al. 2020). This mechanism would explain the inability to change the coat to white in vulgaris least weasels, keeping eumelanin production during the moult. We found no additional associations, and our work could not exclude that winter color and the dorsoventral line patterns are controlled by a single locus, as suggested by Frank (1985). Alternative binding affinities to ASIP isoforms expressed during embryogenesis, which determine melanocyte differentiation and dorsoventral color differences (Millar et al. 1995; Manceau et al. 2011), could also affect the determination of the dorsoventral line. A detailed characterization of the genetic variation along the association region, together with functional assays, will be required to validate this amino acid change as the causal mutation, but our results strongly encourage research in that direction.

A Hard Selective Sweep Underlies the Swedish but Not the Polish Morph Transition

The correlation of the geographic distribution of winter coat color with snow cover variables strongly suggests local adaptation for camouflage (Mills et al. 2018). We thus tested for signatures of selective sweeps, indicative of recent positive selection, along the MC1R genomic region. This region is characterized by significantly reduced nucleotide diversity (in both transition zones) and Tajima’s D (in the Swedish transition only) in vulgaris specimens, compared with the genome-wide vulgaris and nivalis pattern (P < 0.001, 10,000 random samples; fig. 3). Selective sweeps were inferred from the whole-genome sequencing data both using a composite likelihood ratio (CLR) approach derived from genotype likelihoods (Pavlidis et al. 2013) and applying a Hidden-Markov Model approach based on pooled sequencing data (Boitard et al. 2012). The results, controlled for the demography of the populations (inferred from whole mtDNA variation; supplementary table S9, Supplementary Material online), are consistent and support a selective sweep overlapping the association region in vulgaris from Sweden (fig. 3, G, and I) but not in either coloration morph in Poland (fig. 3, H, and J). These results are indicative of a hard selective sweep in Swedish vulgaris but not elsewhere, which agrees with the longer block of linkage disequilibrium found there (fig. 2D). Our genotyping, however, detected two haplotypic blocks in the vulgaris specimens from Poland, both sharing the candidate SNPs but generally differing elsewhere (fig. 2 and supplementary table S5, Supplementary Material online). One of these blocks matches the one identified in Swedish vulgaris, whereas the other occurs predominantly in specimens from the Białowieża Forest (supplementary tables S1 and S5, Supplementary Material online), where the two coloration morphs coexist (Atmeh et al. 2018). We cannot thus discard that a signal of hard sweep could have been eroded or masked in the Polish transition, either by gene flow or by selection acting on multiple variants (soft sweep) (Alachiotis and Pavlidis 2018; Harris and DeGiorgio 2020). Indeed, a hard selective sweep signal overlapping the association region of scaffold 337 was detected when restricting the analysis to the specimens carrying the Białowieża vulgaris block, more predominant in our genomic data set (N = 13; supplementary fig. S8, Supplementary Material online), suggesting that selection may have acted on the vulgaris brown variants in Poland as well. Overall, our results show that past local adaptive pressures favored the vulgaris winter-brown phenotype in Sweden, and possibly in Poland, which is consistent with results obtained in snowshoe and mountain hares (Jones et al. 2018; Giska et al. 2019; Jones et al. 2020b). Considering the Swedish vulgaris specimens, we estimated a rather weak selective sweep (s = 2.1 × 10−5; supplementary fig. S9A, Supplementary Material online), and a similar result was obtained for the partial sweep detected in Polish vulgaris specimens carrying the Białowieża haplotype (s = 1.8 × 10−4; supplementary fig. S9B, Supplementary Material online). The strength of the sweep signals could, however, have been reduced by gene flow, in particular given the dominant nature of the vulgaris variant, which underlies the morph even in heterozygosity (Frank 1985; fig. 2). Given the known impact of the coloration morph in camouflage on environments with different snow covers (Atmeh et al. 2018), this selective sweep may reflect the northward retraction of snow after the last glacial maximum, favoring the spread of the vulgaris winter-brown variant, as previously proposed for snowshoe hares (Jones et al. 2018, 2020a, b). Considering the reported impact of color mismatches in the fitness of seasonal coat color changing species (Imperio et al. 2013; Sultaire et al. 2016; Zimova et al. 2016; Atmeh et al. 2018; Wilson et al. 2018; Melin et al. 2020) and the predicted snow cover decreases under climate change scenarios (Brown and Mote 2009), positive selection on winter-brown alleles may be intensified under future warming, maintaining camouflage under decreased snow cover (Zimova et al. 2016; Jones et al. 2018, 2020a). If looking into the past allows predicting the future, our results underscore the adaptive potential of this standing trait variation to promote rapid evolutionary changes in response to climate change (Mills et al. 2018), even over short temporal frames (Giska et al. 2019).

Fig. 3.

Evidence of a hard sweep on the Swedish but not the Polish population. (A, B) Nucleotide diversity (π) estimated along the first 0.75 Mb of scaffold 337 in nonoverlapping 25 kb windows for each coloration morph from (A) Sweden and (B) Poland. (C, D) Tajima’s D estimated along the first 0.75 Mb of scaffold 337 in nonoverlapping 25 kb windows for each coloration morph from (C) Sweden and (D) Poland. (E–H) SweeD composite likelihood ratio (CLR) estimations along the first 0.75 Mb of scaffold 337 for (E) nivalis and (G) vulgaris specimens from Sweden, and (F) nivalis and (H) vulgaris specimens from Poland. The red dashed line shows the P = 0.01 threshold based on the inferred population history and the gray dashed line is the 99th percentile of the empirical distribution of CLR values. (I) Pool-hmm inference of selective sweeps along the first 0.75 Mb of scaffold 337 in vulgaris specimens from Sweden. Dashed gray lines delimit the association region. (J) Pool-hmm inference of selective sweeps along the entire scaffold 337 in vulgaris specimens from Poland did not identify any putatively selected region. The gray-shaded area in panels (A–H) represents the association region identified from the whole-genome scans.

Conclusions and Broader Implications

Our museomics study linked MC1R to the evolution of adaptive seasonal coat color variation in the least weasel (Mustela nivalis), reinforcing the prominent role of the MC1R–ASIP complex in determining fixed but also seasonal color variation (Våge et al. 2005; Jones et al. 2018; Giska et al. 2019). Even though we could not quantify the expression of the gene in the sampled skins, a limitation resulting from the traditional preservation of biological tissues at museums, our results suggest that the genetic control of seasonal coloration via MC1R in least weasels proceeds from protein-coding variation, which contrasts with the regulatory changes at ASIP influencing the variation of the trait in hares (Jones et al. 2018; Giska et al. 2019). Although the array of taxa with unknown determination of seasonal color variation remains large, these results suggest that the evolution of coloration polymorphisms in seasonal coat color changing species has targets of large effect. However, at the fine scale, the means by which these changes evolved reflect diverse evolutionary solutions that facilitated adaptation to changing climates. The drawers of Natural History Museums are filled with representatives of local adaptation, including variation driving cryptic coloration, and hold an immense potential to unlock a broader understanding of the evolutionary mechanisms driving repeated trait evolution in natural populations. In particular, these resources can be used to identify genetic variation that can anchor adaptive responses to climate change.

Materials and Methods

De Novo Genome Assembly

One road-killed male least weasel (M. nivalis) of nivalis coloration morph was opportunistically sampled in Poland (Grzymały-Sierzputy Stare, along the main road Ostrołęka—Łomża), immediately frozen and stored at −20 °C in RNALater until DNA extraction. Extraction of high-molecular weight DNA from heart tissue and 10× Chromium library preparation were performed at the Genomic Services Laboratory of Hudson Alpha Institute for Biotechnology. Libraries were sequenced in two lanes of an Illumina HiSeq X Ten platform, using a 2 × 150 bp paired-end approach.

The linked-read data were used for de novo genome assembly using the Supernova assembler v.2.1.1 (Weisenfeld et al. 2017). Assemblies with several maximum read number (–maxreads; 1, 1.3, 1.4, and 1.7 billion) were tested, and that with higher N50 was retained (1.7 billion reads). Quality control metrics for the retained assembly were obtained with QUAST v.5.02 (Mikheenko et al. 2018). The completeness of the assembled genome was assessed using BUSCO v.3.1.0 (Simão et al. 2015), based on the mammalia_odb9 database (supplementary table S3, Supplementary Material online).

Sampling and DNA Extraction

In total, 175 M. nivalis specimens were included in this work, originally collected between 1954 and 2010 (supplementary tables S1 and S4, Supplementary Material online). This set aimed at reducing dubious scoring of the nivalis and vulgaris coloration morphs (type I and type II, respectively; Frank 1985; King and Powell 2007) and minimizing putative errors in data records, considering the available material. For least weasels from the collection of the Swedish Museum of Natural History (NRM, Stockholm, Sweden), we examined prepared skins and classified as 1) vulgaris (N = 36 retained for analyses) only specimens with data record indicating original collection between November and March, no evidence of a developing white dorsal coat and a ragged dorsoventral line, and as 2) nivalis (N = 45 retained) specimens with a white dorsal pelage (Frank 1985; King and Powell 2007). Specimens not following these criteria were excluded from the dataset. From Poland, no whole skins were available, and we relied on the original data records of the Mammal Research Institute of the Polish Academy of Sciences (MRI, Białowieża, Poland). We retained 59 specimens, which included 1) 12 classified as nivalis for showing a winter-white coat, and 2) 47 collected in other seasons (between April and October) that were originally classified as nivalis (N = 22) or vulgaris (N = 25) considering the straight or ragged dorsoventral line, respectively (Frank 1985; King and Powell 2007). The dorsoventral line shape has been shown to cosegregate with winter coloration and to possibly be determined by one locus or linked loci (Frank 1985). For the Polish samples, we verified that morph classification respected the known geographic distribution (Atmeh et al. 2018). Finally, 35 specimens from the intermorph crossing experiments conducted by Frank (1985), deposited at the Zoological Research Museum Alexander Koenig (ZFMK, Bonn, Germany), were sampled after examining whole prepared skins. Frank (1985) crossed a nivalis female imported from Sweden first with a vulgaris male from Germany and then with male F1 and first backcrosses (BC1). The original data records and publication (Frank 1985) classified the nivalis/vulgaris coloration morph and confirmed the ability of winter whitening of nivalis specimens but not of the vulgaris individuals. Samples consisted of patches of dry skin. DNA extraction followed Dabney et al. (2013).

Library Preparation and Sequencing

For 79 specimens, genomic DNA libraries were prepared and the whole-genome sequenced: 20 of each color morph from Sweden, and 20 nivalis and 19 vulgaris from Poland. The winter specimens showed at least 95% of brown (vulgaris; and ragged line) or white (nivalis) dorsum (following Mills et al. 2013), and preferentially included those originally collected between December and February, at the height of winter color (King and Powell 2007). For the 36 nonwinter specimens included from Poland, the dorsoventral line classified the morphs (see above). Double-indexed genomic libraries were prepared following Meyer and Kircher (2010) with slight modifications, including the USER enzyme (New England BioLabs) in the blunt-end repair (Kircher et al. 2012) (see supplementary note 1, Supplementary Material online). Indexed libraries were quantified with quantitative polymerase chain reaction (qPCR) and sequenced at low individual coverage on an Illumina HiSeq 1500 (1 × 100 bp single-end) at CIBIO-InBIO’s New-Gen sequencing platform (University of Porto, Portugal) and on an Illumina HiSeq X Ten (2 × 150 bp paired-end) at Novogene Co. Ltd. (Yuen Long, NT, Hong Kong).

Data Processing

Raw sequencing data quality was assessed with FastQC v.0.11.7 (Andrews 2017), and adapter sequences removed with Trimmomatic v.0.38 (Bolger et al. 2014), with filters: TRAILING:15, SLIDINGWINDOW:5:15, MINLEN:25. For paired-end reads, overlapping pairs were merged with PEAR v.0.9.1 (Zhang et al. 2014). Cleaned reads were mapped onto the least weasel de novo genome assembly (scaffolds >50 kb) using BWA-MEM v.0.7.17 with default parameters (Li and Durbin 2009). PCR duplicates were removed with MarkDuplicates from PicardTools v.2.18.17 ( http://broadinstitute.github.io/picard/, last accessed November 16, 2019), and reads were locally realigned with RealignerTargetCreator and IndelRealigner from GATK v.3.8.1 (McKenna et al. 2010). Mapping and coverage statistics were generated with Qualimap v.2.2.1 (Okonechnikov et al. 2016). Specimens with mean genome-wide coverage >2× were randomly subsampled to a maximum of 1.75× to avoid overrepresentation.

To account for potential DNA-sequencing errors from DNA damage of the historical samples (Briggs et al. 2007), in addition to the USER enzyme treatment during library preparation, in silico filtering of C-to-T and G-to-A SNPs was done (Axelsson et al. 2008). Allele frequencies were estimated with ANGSD v.0.923 (Korneliussen et al. 2014), using the following filters (applied in all subsequent ANGSD inferences): -uniqueOnly -remove_bads 1 -minQ 20 -minMapQ 20 -SNP_pval 1e−6. Sites with minor allele frequency ≥0.025 (i.e., 4 out of 158 chromosomes) for T or A were removed from the data set. In addition, only sites with sequencing data for a minimum of six individuals per coloration morph and transition zone were retained for analyses, to decrease biases associated with specimen representation.

Population Structure and Evolutionary History

Population structure was first assessed using a PCA in ANGSD, with the single-read sampling approach (-doIBS 2 -doCov1), appropriate for low-coverage data. SNPs were inferred based on genotype likelihoods (-GL 1 -doMajorMinor 1 -doMaf 1), with the following filters: -minInd 60 -setMaxDepth 285 -doCounts 1 -minMaf 0.025. To minimize nonindependence due to linkage, only SNPs distancing at least 20 kb were considered (filtered using PLINK v.1.90 --bp-space command; Purcell et al. 2007). The proportion of variance explained by each component was calculated with R function eigen. A clustering analysis was then done with NGSadmix v.32 (Skotte et al. 2013), using genotype likelihoods estimated with ANGSD (-GL 1 -doMajorMinor 1 -doMaf 1 -doGlf 2) and setting K = 2. A total of 200 independent runs were optimized in a maximum of 2,000 iterations, and the run with the highest likelihood was retained. A neighbor-joining (NJ) tree was constructed with the BioNJ tree building algorithm of FastME v.2.1.5 (Lefort et al. 2015), based on individual pairwise genetic distances estimated with ngsDist v.1.0.9 (Vieira et al. 2016), from ANGSD genotype posterior probabilities (-GL 1 -doMajorMinor 1 -doMaf 1 -doGeno 8 -doPost 1). TreeMix v.1.13 (Pickrell and Pritchard 2012) was further used to assess the relation of population structure and the nivalis/vulgaris coloration morphs. Allele frequencies for each morph–geography combination (Swedish nivalis, Swedish vulgaris, Polish nivalis, and Polish vulgaris) were calculated with ANGSD (with maximum depth three times the average coverage), using the ferret (M. putorius furo) genome to determine the ancestral state (-GL 1 -doMaf 1 -doMajorMinor 5; see supplementary note 2, Supplementary Material online). The TreeMix input was generated with custom scripts, sampling sites distancing at least 20 kb. The best tree topology was inferred using 500 independent runs, with sample size correction turned-off (-noss option) and a final global round of rearrangements (-global).

Local population structure in the Swedish and Polish least weasel samples was further investigated using independent PCAs. To account for the reduction of sample sizes, ANGSD parameters were adjusted: -minInd 30 -minMaf 0.05 and -setMaxDepth 127 or 158. The number of principal components summarizing meaningful genetic variation (i.e., potentially nonrandom) was inferred via 1,000 randomizations of the observed eigenvalues (Peres-Neto et al. 2005), as implemented in the PCDimension R package (Wang et al. 2018), establishing the significance threshold at P < 0.05. In addition, a pairwise genetic distance matrix was estimated with ANGSD (-doIBS 2 -makeMatrix 1) and significant differences in genetic variation between the coloration morphs were tested using a permutational multivariate analysis of variance (PERMANOVA) on the genetic distance matrix, with the adonis function of the vegan R package (Oksanen et al. 2019).

Mitochondrial DNA Assembly and Phylogeny

Mitochondrial genomes (mtDNA) were reconstructed from the short-read sequencing data, using the de novo assembly approach implemented in NOVOplasty v.3.8.2 (Dierckxsens et al. 2017). Two runs were done per sample, using a kmer size of 21 (paired-end data) or 16 (single-end data), and the sequences of NADH dehydrogenase subunit 1 (ND1) or cytochrome B (CytB) gene as seed, recovered from an available mtDNA sequence of M. nivalis (GenBank accession: NC_020639.1). This sequence was also used as reference to improve assembling of repetitive regions. For specimens with paired-end data, the average insert size estimated from the data was set. Final sequences were determined as the consensus of the independent runs, removing the hypervariable control region, which could be erroneously assembled from short reads. Only sequences longer than 15 kb were retained for analyses (31 sequences from Swedish and 33 from Polish specimens; supplementary table S2, Supplementary Material online).

The mtDNA phylogeny was inferred with BEAST2 v.2.5.2 (Bouckaert et al. 2014), using the MultiTypeTree package (Vaughan et al. 2014) to account for the population structure of the two sampled geographic regions, setting a GTR + I + G site substitution model (selected with jModelTest v.2.1.10 [Darriba et al. 2012], using the Akaike Information Criterion) and a strict molecular clock. Three independent runs of 100 million generations were performed for each data set, sampling trees and parameters every 10,000 generations. Convergence was assessed with Tracer v.1.7.1 (Rambaut et al. 2018), and the final phylogeny was reconstructed based on the combined results of the three replicates, discarding the first 10% as burn-in.

Mitochondrial DNA Demographic Inferences

Given that our individual low-coverage whole-genome sequencing approach did not allow reliably modeling population demography due to expected biases in the discovery of low-frequency SNPs (Buerkle and Gompert 2013), the recovered mtDNA data were additionally used to estimate population history parameters. Effective population sizes (Ne) were estimated with BEAST2 (Bouckaert et al. 2014), using the MultiTypeTree package (Vaughan et al. 2014). Three populations were defined (Sweden, Poland group I, and Poland group II), considering the genetic structure of the Polish samples identified with the phylogenetic inferences (supplementary fig. S3, Supplementary Material online), to prevent biases in Ne estimations. The mutation rate was calibrated using the average corrected divergence (estimated with PAUP v.4.0; Swofford 2003) between the M. nivalis and the M. erminea (GenBank accession: NC_025516.1) mtDNA haplotypes, excluding the control region, and assuming the species split at 3.5 Ma (Koepfli et al. 2008) and a generation time of 1 year (King and Powell 2007; Zub et al. 2009). Population sizes were inferred based on the combination of three independent runs of 100 million generations, discarding the first 10% as burn-in. Further, population size changes in Sweden and Poland were inferred using the coalescent models implemented in BEAST v.1.10.4 (Drummond and Rambaut 2007). For each of the three populations, demographic parameters were estimated under a constant size or an exponential growth tree prior, assuming a strict molecular clock and using the best-fit nucleotide substitution models for that data set. Three independent chains of 100 million generations were run and combined with 10% burn-in. Marginal likelihood estimation (MLE) was conducted under a path sampling/stepping stone sampling (PS/SS) approach for model choice and determined that a constant population size provided the best fit in both sampled geographic regions (log MLE values difference <0.5; Kass and Raftery 1995).

Whole-Genome Scans of Differentiation and Association

Candidate genomic regions to determine the nivalis and vulgaris coloration morphs were first identified using whole-genome scans of differentiation implemented in PoPoolation2 v.1.201 (Kofler, Pandey, et al. 2011). Concatenated bam files per morph were converted to mpileup files with SAMtools v.1.9 (Li et al. 2009), requiring a minimum mapping quality of 20. Resulting files were converted to PoPoolation2 sync format, retaining sites with a minimum base quality of 20 and masking windows of 5 bp around inferred indels. To take advantage of the replicated sampling of color morph transition zones and remove possible confounding effects arising from putative local population structure, a Cochran–Mantel–Haenszel (CMH) test was conducted, identifying sites with replicated allele frequency differences across the independent transition zones. Fisher exact tests (FETs) of allele frequency differences between morphs were also performed per site, in the Swedish and Polish data sets separately and pooled together. Significance thresholds for CMH and FETs were established using the Bonferroni correction for multiple tests. The fixation index (FST) between morphs was also calculated for the Swedish and Polish populations independently, across nonoverlapping sliding windows of 25 kb (with --karlsson-fst correction). Windows with <50% of sites covered and <10 SNPs were discarded. Outlier regions of differentiation were defined based on the 99.9th percentile of FST values when the average window FST ≥ 0.40. Results were robust to other tested window lengths and steps. In all PoPoolation2 runs the following filters were applied: --min-count 4 --min-coverage 10 and --pool-size 40:40 or 40:38 for sliding-window analyses. Maximum coverage per morph was established at two times the mean coverage of each group.

Association tests were performed using genotype likelihoods and case–control scans (-doAsso 1) implemented in ANGSD (Korneliussen et al. 2014), with the same data sets used for FETs estimations. Allele frequencies were estimated based on genotype likelihoods (-GL 1 -doMaf 1 -doMajorMinor 1), and differences summarized with a likelihood-ratio test (LRT). For each contrast, we retained sites represented by ≥50% of the individuals (-minInd 20 or 40) and with maximum depth of two times the mean coverage. An association score test (-doAsso 2) was also run, using the first principal component (PC1) of the global PCA (supplementary fig. S2A, Supplementary Material online) as covariate, to account for possible spurious effects of population structure. The test was run following a dominant model, based on posterior genotype probabilities (-doPost 1) estimated from genotype likelihoods, with the additional filters -minHigh 3 and -minCount 20. Associations were summarized with an LRT, converted to P values assuming a chi-square distribution with one degree of freedom. Significance thresholds were established using the Bonferroni correction for multiple tests.

Regions of elevated differentiation and/or significant association consistently recovered across scans were inspected for gene content. Regions of interest were aligned against the ferret reference genome (Peng et al. 2014) using the blastn application of Blast+ v.2.6.0 (Camacho et al. 2009). Ensembl annotations (https://www.ensembl.org/, last accessed March 13, 2020) and gene predictions from the NCBI database (https://www.ncbi.nlm.nih.gov/, last accessed March 13, 2020) for the ferret genome were used to characterize the genes in these regions.

SNP Genotyping

To validate the association identified in the genome scans, 35 specimens from the original intermorph crossing experiments of Frank (1985), deposited at the Zoological Research Museum Alexander Koenig (ZFMK, Bonn, Germany), were genotyped for selected SNPs. For the candidate region in scaffold 337, 22 SNPs showing high allele frequency differences between morphs along the longer association region found in the Swedish population were selected. In addition, 18 SNPs selected from other windows in the top 0.01% of FST values in the Swedish population (scaffolds 333 and 338) were used as controls. Knowing that coloration morphs segregate in least weasels as nivalis recessive (Frank 1985), selected SNPs had minor allele frequencies <0.1 in nivalis and >0.5 frequency difference to vulgaris, estimated with ANGSD. In a few cases, we relaxed the criteria to consider minor allele frequencies of <0.3 in nivalis and frequency difference >0.3 between morphs. In addition, to analyze segregation patterns of genetic variation across the association in wild specimens and narrow down the associated block, the same SNP set was genotyped in 80 individuals from the Swedish and 58 from the Polish populations. Genotyping was performed using the MassARRAY technology at the Genome Transcriptome Facility in the University of Bordeaux, France. The inferred candidate polymorphisms at positions 165,924 and 165,925 in scaffold 337 were additionally genotyped independently by PCR amplification and Sanger sequencing (see supplementary note 3, Supplementary Material online). The significance of the genotype–phenotype associations was calculated with the SNPassoc R package (González et al. 2007), using the WGassociation function and considering dominant inheritance of the vulgaris allele. The significance threshold of P < 0.05 was corrected for multiple tests using the Bonferroni.sig function.

Sequence Conservation at Candidate Mutation

To assess sequence conservation at the identified missense mutation in MC1R gene, MC1R protein sequences from 100 mammal species (retrieved from NCBI, https://www.ncbi.nlm.nih.gov/, and UniProtKB, https://www.uniprot.org/, last accessed April 21, 2020) were aligned using ClustalW, as implemented in BioEdit v.7.0.5.3 (Hall 1999). The functional impact of the amino acid substitution was then predicted using SIFT (Sim et al. 2012), which takes into account both the conservation of the residue across homologous proteins and the biochemical properties of the amino acid substitution. Functional changes were predicted if the normalized probability of tolerated change was <0.05.

Genetic Diversity and Selection Analyses

Genome-wide nucleotide diversity (π) and Waterson’s theta (θ) were estimated with PoPoolation v.1.2.2 (Kofler, Orozco-terWengel, et al. 2011) for the four coloration morph-geography groups (Swedish nivalis and vulgaris, Polish nivalis and vulgaris). Mpileup files were again produced with SAMtools, with the same filters (see Whole-genome scans of differentiation and association). Diversity indices were estimated in non-overlapping 25 kb sliding windows, using the following filters: --pool-size 40 or 38 --min-count 4 --min-coverage 10 --min-covered-fraction 0.50. Maximum coverage was established at two times the mean of each group. Tajima’s D estimates were also conducted for each color morph group, adjusting the --min-count filter to 2.

Signatures of hard selective sweeps along scaffold 337 were identified using the composite likelihood ratio approach implemented in SweeD v.3.2.1 (Pavlidis et al. 2013). Independent runs were performed for each of the four morph-geography groups, as well as the two haplotype groups identified in vulgaris specimens from Poland (fig. 2). Allele frequencies were estimated with ANGSD from genotype likelihoods, using the filters applied for the TreeMix analysis. SNPs were polarized using the ferret as outgroup (-doMaf 1 -GL 1 -doMajorMinor 5; see supplementary note 2, Supplementary Material online). The ANGSD output was converted for SweeD input format using a custom-made R script, and SweeD was run allowing the inclusion of fixed derived sites and estimating the composite likelihood ratio (CLR) at approximately each 10 kb along the scaffold (--grid 4139). The significance P-value threshold for CLR was defined from distributions generated from 1,000 replicate data sets of 1 Mb for each population simulated using msms v.3.2rc (Ewing and Hermisson 2010), according to the demographic parameters inferred from mtDNA. For the Polish least weasel population, the lower Ne estimate of group I was considered, which provided a conservative threshold. CLR values for each replicate were estimated with SweeD at each 10 kb and P < 0.01 was defined as the CLR significance threshold.

Selective sweeps were additionally inferred with the hidden Markov model approach implemented in Pool-hmm v.1.4.4 (Boitard et al. 2013) for pooled sequencing data (Boitard et al. 2012). The following settings were common to all analyses: -q 20 -e sanger --pred -k 1e-20. For each morph-geography group, minimum coverage was set to 10, whereas for the putative Polish vulgaris haplotypes, half its mean value was used. Maximum coverage was established at two times the mean value, and mean genome-wide theta was set based on PoPoolation estimations. The transition probability between hidden states (k) was determined based on neutral data sets simulated according to the inferred demographic parameters. We simulated 1,000 replicate regions of 1 Mb for each population using fastsimcoal v.2.6.0.3 (Excoffier et al. 2013). Selective sweeps were inferred with Pool-hmm for each replicate data set, using k values ranging from 1e−10 to 1e−30. A conservative value of k = 1e−20, which resulted in a false discovery rate = 0, was retained.

To infer the selection coefficients of the selective sweeps identified in Swedish vulgaris and Polish vulgaris with the Białowieża haplotype, we first ran SweeD on 1,000 bootstrap replicates of scaffold 337, by randomly sampling 50% of the SNPs in the scaffold for each data set. For each bootstrap replicate, the mean alpha (α) value was estimated across a 204 kb region including MC1R and flanking 100 kb regions (coordinates: 337: 64,000–268,000) or a 90 kb region overlapping the Polish sweep signal (coordinates: 337: 220,000–310,000) in each data set respectively, generating a predictive distribution of α. A posterior probability distribution of the selection coefficients (s) was then produced assuming s = r × ln(2Ne)/α (Durrett and Schweinsberg 2005), where r is the recombination rate per base per generation and Ne is the effective population size. Recombination rate was randomly sampled from a uniform distribution based on mouse and human means (0.56–1.26 cM/Mb; Jensen-Seaman et al. 2004), Ne sampled from a uniform distribution based on the limits of the 95% highest posterior density interval from the mtDNA-based demographic inference for the Swedish population or Polish group I (supplementary table S9, Supplementary Material online), and α values sampled from the estimated distribution.

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Click here for additional data file.