Phylogeography of the ermine Mustela erminea and the least weasel M. nivalis from Palae-arctic and Nearctic regions were investigated based on mitochondrial DNA control region sequences. Mustela erminea exhibited a very low level of genetic variation, and geographic structures among populations were unclear. This may indicate that M. erminea recently reoccupied a wide territory in Eurasia following the last glacial retreat. In comparison with M. erminea, genetic variations within and among populations of M. nivalis were much greater. Molecular phylogenetic relationships showed that two lineages of M. nivalis occurred in the Holarctic region: one spread from the Eurasian region to North America, and the other occurred in south-eastern Europe, the Caucasus and Central Asia. The results suggest either mitochondrial DNA introgression among populations of south-eastern Europe, the Caucasus and Central Asia, or ancestral polymorphisms remaining in those populations. Contrastive phylogeographic patterns between the two mustelid species could reflect differences of their migration histories in Eurasia after the last glacial age.
The ermine Mustela erminea Linnaeus, 1758, inhabits the tundra, boreal forest and deciduous forest zones of the Palaearctic and Nearctic regions: Europe including Britain and Ireland but excluding most parts of Mediterranean zones; Siberia and Asia including the Caucasus, Himalayas, Mongolia, Manchuria, and Far East including Sakhalin, Hokkaido and Honshu of Japan. This species is also distributed across North America, south to California in the west and to Maryland in the east (Corbet, 1978; Nowak, 1991; Wozencraft, 1993). Although many subspecies of M. erminea have been recognized in Eurasia and North America (Ellerman and Morrison-Scott, 1951; Hall, 1951; Heptner et al., 1967; King, 1983), an unequivocal taxonomic classification has not yet been made. Corbet (1978) recognized M. erminea throughout Eurasia and North America as one subspecies, and only insular populations from Ireland and the Honshu Island of Japan as separate subspecies.
The least weasel Mustela nivalis Linnaeus, 1766, also has a circumboreal range throughout the Holarctic region, including Europe, North Africa, and a major part of Asia and North America (Corbet, 1978; Nowak, 1991). Over this wide area, M. nivalis displays such a high degree of morphological variations that the systematic relationships among different geographic forms are unresolved. The most recent revisions of the morphometric variation of M. nivalis were carried out by van Zyll de Jong (1992), Reig (1997) and Abramov and Baryshikov (2000).
Thus, Mustela erminea and M. nivalis are largely sympatrically distributed over the major parts of Eurasia and North America. Paleontological evidences indicate that the two species have been typical small predators of the same biotic association since the Middle Pleistocene (Kurten, 1968; Kurten and Andersen, 1980). The two species are morphologically and ecologically very similar to each other. Therefore, both species were believed to be closely related and they were formerly placed in one subgenus Mustela (Ellerman and Morrison-Scott, 1951; Youngman, 1982; Nowak, 1991; Wozencraft, 1993). Their common circumpolar distribution suggests that they might have experienced the similar history in the post glacial recolonization. However, comparative analyses of geographical variations in morphological characters show differences between the two species: the geographic variability of morphology in M. nivalis is larger than that of M. erminea (Heptner et al., 1967; Ralls and Harvey, 1985; Meia 1990; Meia and Mermod, 1992).
Recently, Kurose et al. (1999) found nucleotide sequence differences in the mitochondorial DNA (mtDNA) control region between the Hokkaido and Honshu populations of M. erminea and between those of M. nivalis in Japan. The mtDNA (Kurose et al., 1999) and karyotypic (Obara, 1991) differentiations of M. nivalis could have been fixed to the Honshu Island population, after their geographic isolation in the Middle Pleistocene. Meanwhile, the mtDNA control region sequences did not show any clear divergence between the Hokkaido and Honshu populations of M. erminea (Kurose et al., 1999). Fleming and Cook (2002) examined DNA sequence variations in the mitochondrial cytochrome b gene of M. erminea in Alaska and found three lineages. One lineage includes southeast Alaskan endemics and ermines from western Canada and the coterminous United States. The second lineage appears restricted to the Prince of Wales Island complex in southeast Alaska, the Graham Island and British Columbia. The third lineage has a Beringian distribution including Alaska and north-eastern Russia (Magadan).
In the present study, we compared phylogeographic patterns of mtDNA control region sequences between the two sympatric mustelid species across Eurasia and North America, and then discuss their post glacial history and the subspeciation.
MATERIALS AND METHODS
Animals and DNA extraction
Profiles and sampling localities of animals examined in the present study are described in Table 1 as well as Figs. 2a and 3a. All samples were muscle tissues preserved in 70% ethanol for approximately 10–160 years. Total DNAs were extracted from homogenates of tissues (approximately 1×1×1 cm), according to the phenol/proteinase K/sodium dodecyl sulphate method of Sambrook et al. (1989) with some modifications (Masuda and Yoshida, 1994). DNA extracts were concentrated using Centricon-30 microconcentrators (Amicon), because DNAs contained in long-preserved tissues were fragmented. The STE buffer (100mM NaCl/10mM Tris pH 7.5/1mM EDTA) treated by the same procedure was used as a negative control in subsequent polymerase chain reaction (PCR) amplification. MtDNA control region sequences (Accession numbers AB006717, AB006721, AB006722, AB006727–AB006733) from Japanese populations of M. nivalis and M. erminea (Kurose et al., 1999) were added to the present analysis. The sequence (AB007327) of Mustela itatsi was used as an outgroup, because this species is phylogenetically closer to M. erminea and M. nivalis in the genus Mustela (Kurose et al., 2000).
Profiles of ermines Mustela erminea (MER) and least weasels Mustela nivalis (MNI) examined in the present study
PCR, sequencing and data analysis
Since DNAs in long-preserved tissues were fragmented, four portions of the mtDNA control region were amplified separately using two primers from those described in Table 2 and Fig. 1, except for MNE-D8 and MNE-D9 which were used for the following sequencing. From sequences reported by Shield and Kocher (1991) and Árnason and Johnsson (1992), primer Cb-z was modified and primer D4 was designed, respectively (Kurose et al., 1999). Primers DS1, MSD and MER-R were cited from Kurose et al. (1999). The other primers (Tasble 1) were newly designed in the present study based on control region sequences of Japanese populations of M. nivalis and M. erminea, previously reported by Kurose et al. (1999). Combination of the four portions yielded approximately 600 base-pairs (bp) of the 5′ end of the mtDNA control region, which were reported to be polymorphic in Japanese populations of these species (Kurose et al., 1999). Symmetric PCR amplification to generate double-stranded DNAs was done in a total volume of 50 μl of the reaction mixture including 20 mg of bovine serum albumin (Boehringer). PCR included 35 cycles using a DNA thermal cycler (PJ2000, Perkin-Elmer Cetus) as follows: 94°C, 1 min denaturing; 50°C, 1 min annealing; 72°C, 1 min extension, and then the reaction was completed at 72°C for 10 min. To varify positive amplification, an aliquot of 10 μl from the PCR products was electrophoresed on a 2% agarose gel, stained with ethidium bromide, and visualized under an ultraviolet illuminator. The remaining 40 μl of the PCR product was purified with the QIAquick purification kit (Qiagen). Purified PCR products were labelled using the dye-terminator pre-mix kit (Amersham) on a DNA thermal cycler (PCR system 9700, Perkin-Elmer) following the manufacturer's instruction and sequenced using the ABI Prism™ 377 DNA sequencing system (Perkin-Elmer). For sequencing, all PCR primers except for Cb-z and D4 were used.
Primers for amplification and sequencing of the mtDNA control region, which were newly designed in the present study
Sequence alignment was done using GeneWorks (Intelligenetics). Insertions and/or deletions (indels) were compensated through observation by eyes and eliminated for sequence analysis. Neighbor-joining trees (Saitou and Nei, 1987) using MEGA (Kumar et al., 1993), based on genetic distances of Kimura's (1980) two-parameter model, were constructed. The maximum parsimonious method was done using PAUP* 4.0b10 (Swofford, 1998).
Samples with successful sequencing results are listed in Table 1. However, the table does not include other samples, where no PCR products were obtained because of DNA degradation in aged specimens.
Sequence variations in Mustela erminea
An alignment of the mtDNA control region sequences of Mustela erminea (Table 3) showed 40 sites to be variable among 579 bp (excluding one indel at nucleotide site 71) with the transitional bias. Intraspecific sequence diversity (percentage differences) varied from 0.17% to 4.32% (1.16% in average). The neighbor-joining phylogenetic tree (Fig. 2b) clearly showed a low level of genetic differentiation between populations of M. erminea in Eurasia. Most of the clusters of haplotypes were not supported with bootstrap values greater than 71% (Fig. 2b). Two haplotypes E15 (RSEK2) and E5 (RSEK4) from eastern Kazakhstan (locality No. 3 in Fig. 2a) were positioned separately from each other. The two haplotypes E17 (RANA1) and E7 (RANA2) from north-eastern Siberia were also remotely related to each other (Table 1, Fig. 2b and 2c). The haplotype E11 (TOG1 and YAM1) from the Honshu Island of Japan (locality No. 12 in Fig. 2a) was most closely related to the haplotype E10 (RKAM1) from Kamchatka (locality No. 8 in Fig. 2a), although the relationship between them was relatively weak (23% bootstrap value). In addition, the haplotype E18 (US3) from North America (locality No. 15 in Fig. 2a) was closer to E17 (RANA1) from north-eastern Siberia (locality No. 7 in Fig. 2a). Thus, the phylogenetic relationships among haplotypes did not relate closely to their geographic distribution throughout Eurasia. Genetic distances of North American individuals after divergence were much larger than those observed in Eurasian populations (Fig. 2b). The phylogenetic relationships of haplotypes by maximum parsimony (Fig. 2c) were similar to those by the neighbor-joining relationships (Fig. 2b).
Sequence variation of the mtDNA control region (581–582 bp) in the ermine Mustela erminea
Sequence variations in Mustela nivalis
A sequence alignment of the mtDNA control region sequences of Mustela nivalis (Table 4) showed 61 sites to be variable among 579 bp (excluding seven indels) with the transitional bias. Intraspecific sequence differences varied from 0.17% to 6.04% (2.91% in average). The neighbor-joining tree indicated that individuals of M. nivalis were divided into seven clusters A–G with 55–100% bootstrap values (Fig. 3b).
Sequence variation of the mtDNA control region (581–591 bp) in the least weasel Mustela nivalis
Sequence differences in Cluster A with a 55% bootstrap value were between 0.17% to 1.21%. Cluster A includes haplotypes N1 (RCAU1), N2 (RROS1), N4 (RCAC1) and N5 (RETU1) from the Caucasus and haplotypes N3 (RASK1) and N2 (RASK2) from southern Ukraine and haplotypes N1 (RTUR2) from Turkmenistan. Cluster B with a 79% bootstrap value showed sequence differences varying from 0.17% to 0.35%. Three specimens (RALT1, RUMA1, and RLEN3) shared an identical haplotype N7, despite these individuals originating from distant localities (see Fig. 3a). As reported by Kurose et al. (1999), cluster C comprising the Honshu Island population of Japan was defined by one haplotype N9, which was shared by three individuals (AKI1, IWA1, and IWA2), and cluster D (the Hokkaido Island population of Japan; 87% bootstrap value) showed intra-cluster sequence differences varying from 0.35% to 1.04%. Cluster E includes haplotypes N14 (US1) and N15 (US2) from North America with a 100% bootstrap value. Cluster F (one haplotype N16; RKAR2 from Turkmenistan) was split from clusters A–E. Cluster G, with a 100% bootstrap value, showed intra-cluster sequence differences varying from 0.52% to 1.55% and consisted of individuals from distant localities: N17 (RCRY1) from Crimea, N18 (RCAU2) from the Caucasus, N19 (RKAR1) from Turkmenistan and N20 (RJIR6) from western Kazakhstan (Fig. 3b). Within-population nucleotide diversities of M. nivalis were higher than those of M. erminea mentioned above. For instance, the population from eastern Turkmenistan contained very distinct haplotypes N19 (RKAR1), N16 (RKAR2) and N1 (RTUR2) that belonged to three clusters A, F and G, respectively. Two haplotypes N1 (RCAU1) and N18 (RCAU2) of the morphologically similar specimens from the same locality in the Caucasus were included in different clusters, A and G, respectively. The phylogenetic relationships analyzed by maximum parsimony (Fig. 3c) were similar to those analyzed by the neighbor-joining method (Fig. 3b).
Genetic population structure of Mustela erminea
Genetic differentiation among populations of Mustela erminea was relatively lower than that of M. nivalis. Phylogenetic relationships based on sequence differences did not correspond with geographic distances of sampling localities (Fig. 2b and 2c). For instance, Japanese haplotypes were separated into two groups and positioned with clusters of the Eurasian continental haplotypes. North American haplotypes were also clustered into two groups and separately positioned with Eurasian haplotypes.
Phylogenetic data based on cytogenetic (Graphodatsky et al., 1976; Obara, 1991) and mtDNA (Kurose et al., 2000) analyses suggest that M. erminea was first split from other species of the genus Mustela and that this species has a more ancestral form among the genus Mustela. Meanwhile, M. erminea likely appeared in central Europe in the early Middle Pleistocene and it was one of the common small carnivores in central and eastern Europe during the Middle Pleistocene and later (Kurten, 1968). The original European population of M. erminea could have colonized eastern parts of the Palaearctic region and crossed the Bering land-bridge into North America, where the earliest known remains are of the Illinoian age (Kurten and Andersen, 1980). Although recent M. erminea slightly differs from the Early Pleistocene and Late Pliocene ermine-like forms, the geographical variability in cranial and external characters of the Palaearctic and Nearctic ermines is also very limited (Ralls and Harvey, 1985; Eger, 1990; Meia, 1990; Meia and Mermod, 1992). The results of the present study may indicate that the intercontinental genetic difference between Eurasian and North American populations is not large, in congruence with the low level of intra-specific variation obtained in previous morphological studies. One explanation for the low variation in M. erminea is that colonization of Eurasia might have occurred quickly during a short period after the last glacial age. After the original colonization by the ancestral populations of eastern Palaearctic regions including Japan, there may have been continuous exchanges among populations, resulting in homogenization of mtDNA. Meanwhile, it was reported that populations of the gray wolf Canis lupus showed a similar pattern of mtDNA control region sequences, that is, low genetic variations and no geographic structures in Eurasia (Tsuda et al., 1997; Vila et al., 1997). It was suggested that mtDNA genotypes were randomly fixed to be monomorphic because of decreasing and fracturing of available habitats for wolves and a decrease in the size of their populations. Mustela erminea might also have experienced the similar population history.
Meanwhile, the phylogenetic trees indicate that the mutation rates of the North American haplotypes (E8, E9 and E18) are likely more rapid than those of Eurasian haplotypes (Fig. 2b and 2c). North American haplotypes were positioned not outside but within clusterings of the Eurasian haplotypes. The future examination by comprehensive sampling in North America may reveal the genetic specificity of the mtDNA control region of North American M. erminea. The differences found in the present study between North American ermines from Alaska (Kodiak and Kenai) and inner USA (Minnesota) correspond with data of Fleming and Cook (2002). These authors found three distinct lineages in North America – Beringian (includes specimens from northeastern Alaska and north-eastern Russia), the Continental (distributed from the south-eastern Alaska throughout western Canada and across the coterminous USA), and the Island lineage (found on a few islands of southeastern Alaska).
Genetic population structure of Mustela nivalis
In comparison with M. erminea, the genetic variation within and among populations of M. nivalis was much larger. Both neighbor-joining and parsimonious methods produced the similar phylogenetic relationships among haplotypes (Fig. 3b and 3c). There are clearly two distinct mtDNA lineages for M. nivalis: one lineage comprizing clusters A–F with 93–95% bootstrap values (Fig. 3b and 3c) is distributed across Eurasia including Japan and North America, and the other lineage (cluster G) with a 100% bootstrap value (Fig. 3b and 3c) was found in some populations of the Caucasus, Central Asia (Turkmenistan, Kazakhstan) and south-eastern Europe (southern Ukraine).
Recently, Abramov and Baryshnikov (2000) investigated cranial characters, sizes, proportions of the body and tail, and the coloration of pelage, and reconstructed subspecific taxonomy of M. nivalis as well as possible phylogenetic relationships. The results of the present study mostly support their phylogeny.
According to the above-mentioned taxonomy, most parts of northern Eurasia (including the Hokkaido Island of Japan) were occupied by the small-sized subspecies M. n. nivalis. All specimens from clusters B and D (European parts of Russia, Altai, Eastern Siberia, Hokkaido) belong to this form.
Obara (1991) reported the occurrence of the Robertsonian translocation of chromosomes in the Honshu Island population of Japan (2n=38), which distinguishes this population from the other Eurasian (including Hokkaido Island of Japan, 2n=42) and North American populations (2n=42). Based on the karyological difference, Obara (1991) considered the Honshu Island population as a distinct species M. namiyei Kuroda, 1921. However, genetic distances between the two populations estimated from the mtDNA cytochrome b (Masuda and Yoshida, 1994) and control region (Kurose et al., 1999) did not support the specific level of divergence. The molecular phylogeny obtained in the present study again supports that the Honshu population is genetically differentiated but is not more than the subspecific level. Interestingly, the population from Honshu Island (cluster C) was closer to the continental M. n. nivalis (cluster B) (with 55–61% bootstrap value) than the Hokkaido Island group (cluster D) (Fig. 3b and 3c).
Moreover, the Hokkaido Island group (cluster D) was closely related to the North American group (cluster E), while this was supported with <50% bootstrap values (Fig. 3b and 3c). North American specimens used in the present study were classified as M. n. allegheniensis, which is morphologically closer to the nominotypical subspecies M. n. nivalis that also includes weasels from Hokkaido (Abramov and Baryshnikov, 2000). Based on craniological data, Reig (1997) considered the North American least weasel as distinct species M. rixosa (Bangs, 1896). However, the present molecular data (Fig. 3b and 3c) as well as the results of previous morphological analysis (Abramov and Baryshnikov, 2000) do not support the specific level of divergence between Eurasian and North American least weasels.
The other three clusters (A, F and G) consisted of medium-sized weasels from southern Ukraine (subspecies vulgaris) and the Caucasus (subspecies caucasica), and large-sized weasels from Turkmenistan and western Kazakhstan (subspecies heptneri) (subspecific ranks according to Abramov and Baryshnikov, 2000). The clustering of hap-lotypes in the phylogenetic tree (Fig. 3b and 3c) did not correspond with geographically and morphologically expected relationships between populations. Nor did the clustering always relate to body size and color patterns. The least weasels from these localities shared the haplotypes from two different lineages. The results of the present study indicate that mtDNA introgression may have occurred between populations distributed in the Caucasus, southern Ukraine and Central Asia, or polymorphic status of their ancestral populations still remains in those areas. Hewitt (1999) reported that many animal species of Europe show lower genetic diversity in northern populations that have expanded rapidly. By contrast, southern populations in refugial regions show relatively greater diversity.
The ancestry of the recent M. nivalis can be traced back to M. praenivalis of the Late Pliocene and Early Pleistocene in central Europe (Kurten, 1968). Probably more ancestral forms of M. nivalis were large long-tailed weasels of northwestern Africa, southern Spain and Mediterranean (form numidica) (Frank, 1985; Abramov and Baryshnikov, 2000). In the past, these large weasels were more widely spread in Europe and might also penetrate western parts of Asia. Small weasels of the form nivalis appeared later, probably in the boreal regions of the Palaearctic. The sympatric occurrence of fossil remains of large and small weasels in cave layers is known for the Late Pleistocene sites of the Caucasus, Poland, Slovakia and Crimea (see Abramov and Baryshnikov, 2000 and references therein). These small weasels could have penetrated North America through Beringia. Our results suggest that the form nivalis might have colonized the most parts of the Eurasian continent during a short period, probably after the last glacial age. The medium-sized weasel (vulgaris-boccamela) was evidently formed in eastern parts of the Mediterranean region during the Late Pleistocene. From southern refugia, weasels of this form could have recolonized most parts of Europe, displacing nivalis to the north and east (see Kratochvil, 1977). To examine this process, further studies are required by using more specimens from the different areas including the Mediterranean and western Europe.
The present study revealed that the genetic variations between Palaearctic populations were much larger in M. nivalis than in M. erminea. This contrasting result is in congruence with the variability of cranial and external characters of M. nivalis which is much larger than that of M. erminea (Heptner et al., 1967; Ralls and Harvey, 1985; Meia 1990; Meia and Mermod, 1992). These findings reflect that M. nivalis has a wider distribution and inhabits more diverse biotopes and possibly has a higher adaptability to environments than M. erminea.
We are grateful to the Institute of Low Temperature Science of Hokkaido University for technical support. This study was supported in part by Grants-in-Aid for Scientific Research and the 21th Century COE Program “Neo-Science of National History” at Hokkaido University financed from the Ministry of Education, Culture, Sports, Science and Technology, Japan.