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

Table 1.

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.

Table 2.

Primers for amplification and sequencing of the mtDNA control region, which were newly designed in the present study

Fig. 1.

Schematic diagram of the mtDNA control region of Mustela erminea and M. nivalis. Arrows show primers used for PCR amplification and sequencing. Detailed information is described in ‘MATERIALS AND METHODS’ in the text and Table 1.

Fig. 2.

(a) Sampling localities of Mustela erminea in Eurasia. The small map at the upper right shows the sampling localities in North America. Locality numbers refer to those in Table 1. (b) Neighbour-joining relationships among haplo-types (581–582 bp) of M. erminea mtDNA control region. The homologous sequence of M. itatsi was used as an outgroup. The scale indicates genetic distance estimated with Kimura's (1980) two-parameter model, excluding indel sites. Numbers (%) on internal branches are bootstrap values derived from 1,000 replications. Haplo-type numbers, locality numbers in brackets and locality names in parentheses refer to those in Table 1 and Fig. 2a. (c) Maximum parsimony tree using haplotypes of M. erminea. The sequence of M. itatsi was used as outgroup. The number ‘53’ on the branch indicates a bootstrap value from 100 replications, however, low bootstrap values less than 50% are not shown. The horizonal branch length is proportional to distances as shown by the scale bar (number of substitutions). Haplotype numbers, locality numbers in brackets and locality names in parentheses refer to those in Table 1 and Fig. 2a.

Fig. 3.

(a) Sampling localities of Mustela nivalis in Eurasia. The small map at the upper right shows the sampling localities in North America. Locality numbers refer to those in Table 1. (b) Neighbour-joining relationships among haplo-types (581–591 bp) of the Mustela nivalis control region sequences. The sequence of M. itatsi was used as an outgroup. The scale indicates genetic distance estimated with Kimura's (1980) two-parameter model, excluding indel sites. Numbers (%) on internal branches are bootstrap values derived from 1,000 replications. Haplotype numbers, locality numbers in brackets and locality names in parentheses refer to those in Table 1 and Fig. 3a. ▴, Cluster A; •, Cluster B; ▾, Cluster C; ♦, Cluster D; #, Cluster E; ★, Cluster F; ▪, Cluster G. (c) Maximum parsimony (MP) tree reconstructed using haplotypes of M. nivalis. The homologous sequence of M. itatsi was used as outgroup. Numbers indicate bootstrap values from 100 replications, however, low bootstrap values less than 50% are not shown. The horizonal branch length is proportional to distances as shown by the scale bar (number of substitutions). Haplotype numbers, locality numbers in brackets and locality names in parentheses refer to those in Table 1 and Fig. 3a.

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

Table 3.

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

Table 4.

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.