The population size of the sika deer Cervus nippon on Hokkaido Island of Japan had been remarkably reduced because of heavy hunting pressure since the beginning of Meiji Period and effects of heavy snow in 1879 and 1881. After that, the number of sika deer in Hokkaido has increased gradually due to the protection by the Hokkaido government. In the present study, in order to investigate the bottleneck effects, we analyzed ancient mitochondrial DNA (mtDNA) on sika deer bones excavated from archaeological sites just before Meiji Period. On 86 of 113 bones from 13 archaeological sites of Ainu Culture Period (17–19th centuries), 602 base-pair fragments of the mtDNA control region were successfully sequenced. Consequently, we found three new haplotypes (g-, h- and i-types) which had not been identified in modern sika deer. In addition, four haplotypes (a-, b-, c- and d-types) identified from modern sika deer were also found in the archaeological deer. The new haplotypes and previously reported hapoltypes from sika deer of Hokkaido were phylogenetically much closer to each other, compared with those of modern sika deer from Honshu, Kyushu and the Chinese continent. Geographical distribution patterns of haplotypes of the ancient population were different from those of the modern population in Hokkaido. Our findings indicated that their genetic diversity was reduced through the bottleneck and that population structures of sika deer were changed widely in Hokkaido due to genetic drift.
Cervus nippon yesoensis, a subspecies of the sika deer C. nippon occurring in Japanese islands and eastern Asia, is a population on Hokkaido Island of Japan, whose body sizes are the largest among the subspecies (Abe et al., 1994). Before developing Hokkaido in the beginning of Meiji Period, a large population of the sika deer occurred throughout Hokkaido. Because of heavy hunting pressure since the beginning of Meiji Period and effects of heavy snow recorded in 1879 and 1881, however, the population size of the sika deer was remarkably reduced in Hokkaido (Inukai, 1952). After that, the number of the sika deer increased gradually due to the protection by the Hokkaido government and currently their overpopulation makes a great damage to the agriculture and forestry (Kaji, 1995; Hokkaido Government, 1986, 1994).
Nagata et al. (1999) reported that C. nippon occurring on Japanese islands is divided into northern and southern groups, based on molecular phylogenetic analysis of the mitochondrial DNA (mtDNA) control region. The population of Hokkaido was included in the northern group, which could have immigrated to Japanese island, through landbridges, from the continent in the late Pleistocene (approximately 20,000–30,000 years ago) (Harunari, 2001). Goodman et al. (2001) studied microsatellite polymorphisms and reported that genetic diversity of the northern group was lower than that of the southern group in Japanese islands. Especially in the Hokkaido population, genetic diversity on mtDNA and nuclear microsatellites is remarkably low, compared with that of sika deer of Honshu and Kyushu, possibly because of bottleneck effects by the heavy snow in Meiji Period (Nagata et al., 1998a, 1998b).
Since DNA data from archaeological remains provide invaluable information on phylogeography and history of Japanese mammals as reported in brown bears (Masuda et al., 2001) and pigs and wild boars (Watanobe et al., 2001), ancient DNA analysis was considered to be useful to study the changing process between ancient and modern deer populations in Hokkaido. If genetic population structures of sika deer before Meiji Period is restored, bottleneck effects on the population due to the heavy snow and hunting pressure will be more clarified.
In the present study, we analyzed ancient mtDNA from sika deer bone remains which were excavated from archaeological sites before Meiji Period in Hokkaido. Then, we compare the ancient genetic features with those of modern sika deer in Hokkaido, and discuss changes of genetic population structures due to bottleneck effects.
MATERIALS AND METHODS
A total of 113 bone remains of sika deer examined in the present study were excavated from archaeological sites in Hokkaido (Table 1, Fig. 1). Deer bones of Bibi were excavated from a trace of a deer meat factory which was closed in spring of 1880: one year before (1879) and next year (1881) there was heavy snow in records (Inukai, 1952). Many of those bones were crushed with probable iron apparatus, indicating that organs including bone marrows were utilized in the factory. The other sites were archaeologically determined to be developed at Ainu Culture Period (17–19th centuries). To avoid analyzing bone remains from identical individuals, we used parts of bones in the same positions and the same sides, or with different sizes, or from different strata from one archaeological site.
Profiles of bone remain samples analyzed in the present study
DNA extraction from bone remains
The DNA extraction from bone remains was carried out according to the method of Masuda et al. (2001). Parts of bones were powdered with an electric drill. Approximately 0.5 g of powders were suspended with 5 ml of 0.5M ethylenediamine tetraacetic acid (EDTA) in a 15 ml-plastic tube with rotating at room temperature. The EDTA solution was changed every hour and centrifuged at 3,000 rpm for 15 minutes until the supernatant became clear. The EDTA solution was usually changed twelve times per sample, and then samples with the last EDTA solution were rotated overnight at room temperature. The last sediments were resuspended in 5 ml of 0.5M EDTA and 100 μl of 10 mg/ml proteinase K, and incubated at 37°C overnight with rotating. The solution was extracted through the phenol-chloroform extraction method (Sambrook et al., 1989). Then, the extracts were concentrated into approximately 50 μl of TE buffer using VivaSpin 6 Concentrator (Sartorius). Through experiments, disposable apparatuses such as plastic tubes, tips, and gloves were used to prevent contamination of external DNA. To monitor possible DNA contamination, a blank extraction without bone powders was performed using the same procedure as mentioned above.
Polymerase chain reaction (PCR) amplification
Because of fragmentation of aged DNA, a 602 base-pair (bp) fragment of the mtDNA control region of ancient sika deer were separated into three portions (Fragments 1-3, Fig. 2), and we amplified them with PCR primers reported by Nagata et al. (1998a) and other primers newly designed in the present study (Table 2). Figure 2 shows positions of primers for PCR and the following direct sequencing. Programs of PCR were as follows: one cycle of 94°C for 3 min, and then 45 cycles of 94°C for 1 min, 50°C for 1 min, and 72°C for 1 min in a total volume 50 μl of a reaction mixture containing 1 μl of DNA extracts, 5 μl of 10×buffer, 4 μl of dNTP, 0.5 μl of rTaq DNA polymerase (Takara), each primer at 0.5 μl, and 20 μg of bovine serum aibumin (BSA, Boehringer). The BSA was used to eliminate some effects of PCR inhibitors which were often contained in aged specimens. One cycle of 72°C for 10 min was performed for the reaction completion. If the first PCR failed in amplification, the secondary PCR, semi-nested PCR or nested PCR was carried out with 30–50 cycles. In this case, BSA was not added in reaction mixtures. The primer pair of LD5 and HD12 was used for PCR amplification of Fragment 1. When PCR failed with LD5/HD12, the primer pair LD5in/HD12 was used. On Fragment 2, the primer pair LD11/HD8 or LD13/HD10 was used for the first PCR. When PCR failed with LD11/HD8, the primer pair LD13/HD10 was used for nested PCR. If the first PCR failed with LD13/HD10, the primer pairs LD13in/HD10 and LD13/HD10in were used for semi-nested PCR. On Fragment 3, the primer pair of LD15 and HD2 was used for PCR amplification. When PCR failed with LD15/HD2, the primer pair LD15/HD2in was used. All PCR amplifications were carried out in a PCR thermal cycler TP400 (Takara). The PCR products were purified using the centrifugal dialysis kit QIAquick (Qiagen), and used as template for the direct sequencing.
Primers for PCR amplification and DNA sequencing of the Hokkaido sika deer mtDNA control region
Direct sequencing of PCR products
Sequencing reaction was done with the Thermo Sequence pre-mixed cycle sequencing kit (Amersham). Sequencing was carried out using an automated sequencer Hitachi SQ-5500.
Genetic distances between haplotypes were calculated by the Kimura's (1980) two-parameter method. Gap sites were eliminated from calculation of genetic distances. Phyogenetic trees were reconstructed using the neighbor-joining method (Saitou and Nei, 1987) in MEGA (Kumar et al., 1993). Bootstrap values (Felsenstein, 1985) were derived from 1,000 replications. Six haplotypes of modern deer of Hokkaido (a- to f-types: accession nos. D50128, D50129, AB004297- AB004300; Nagata et al., 1998a) and haplo-types of modern deer of Honshu, Kyushu, small islands and China (AB012366-AB012384; Nagata et al., 1999) were quoted. The sequence (AB012385) of the red deer Cervus elaphus was used as outgroup.
MtDNA haplotypes from archaeological deer
Fragments 1–3 of the mtDNA control region (Fig. 2) were successfully PCR-amplified from 86 of 113 sika deer bone remains of Hokkaido, and a total of 602 bp was sequenced for each of the 86 samples. Compared with six mtDNA haplotypes (a-, b-, c-, d-, e- and f-types) of modern sika deer of Hokkaido reported by Nagata et al. (1998a), three new haplotypes (named g-, h- and i-types) were identified from three archaeological samples examined in the present study (Table 3). In addition, a-, b-, c- and d-types were also found from 83 of 86 samples. Because all of five bone samples from Aoshimanai Site, where individuality among the samples was not morphologically clear, shared the c-type, we regarded them as one individual with the c-type at this site. In addition, on Furetoi samples, where individuality was unclear, the a-type was identified from two bones and the c-type was from one bone. Therefore, we regarded them as one individual with the a-type and one with the c-type from this site. Consequently, we identified mtDNA haplotypes from 81 samples. The three newly identified sequences will appear in the DDBJ nucleotide database ( http://www.ddbj.nig.ac.jp/Welcome.html) with the following accession numbers: AB118753 (g-type), AB118754 (h-type), AB118755 (i-type). No successful results were obtained from 10 samples at Kushiro-katsurakoi-fushikokotan-chashi Site and 10 samples at Akan-shimoninishibetsu Site (Table 1), because of probable degradation of remained DNA.
Frequencies of mtDNA haplotypes in 11 archeaological sites of Hokkaido
Table 4 shows nucleotide substitutions between haplo-types reported by Nagata et al. (1998a) and those identified in the present study. In modern sika deer of Hokkaido, transitional substitutions were found at four sites (A/G at nt 178, 251 and 499; C/T at nt 457). Three new haplotypes (g-, hand i-type) were identified only from Usu-oyakotsu 1 Site. In the g-type, a new combination of substitutions was seen at the four sites (nt 178, 251, 457, 499), and one transition (T>C) occurred at a new nucleotide site, nt 267. The other sequences in the g-type were the same as a-type. The h-type had two new transitions: C->T at a new site, nt 305, G>A at another new site, nt 376. The other sequences in the h-type were the same as d-type. Meanwhile, the i-type had eight motifs of approximately 40 bp which consisted of a duplicate of the fifth motif in the a-type sequence, while the other ancient deer shared seven motifs in the mtDNA control region (Fig. 2). No transversions were found among all ancient and modern deer.
Nine haplotypes of the mtDNA control region of sika deer in Hokkaido
Distribution patterns of mtDNA haplotypes in Hokkaido
Figures 1a, 1b and 1c show geographical distribution patterns of a- to c-types in modern deer (Nagata et al., 1998a) and ancient deer of the present study. The frequency of the a-type was 58 of 81 samples (71.7%) and highest at most archaeological sites (Fig. 1a, Table 3). Also in the modern deer, the a-type was major mainly around coniferous forests of the Akan area (Nagata et al., 1998a). The b-type frequency from archaeological deer was 3.7% (3/81 animals) and lower than that from modern deer. In the modern population, the b-type was localized in northern parts of the Daisetsu area, while this type was found in ancient deer at different areas such as Bibi Site and Benten Site in Ishikari-lowlands (Fig. 1b). Meanwhile, the c-type which occurs mainly around the Hidaka area in the modern population was found in ancient deer of different areas: Bibi Site and Benten Site in Ishikari-lowlands and Onnebetsu Site in eastern Hokkaido (Fig. 1c). Three deer from Benten Site, Penakori 1 Site and Opaushinai 1 Site shared the d-type. The e- and f-types were not identified from ancient deer, while they were found from modern deer.
Phylogenetic relationships among mtDNA haplotypes
Fig. 3 indicates most parsimonious relationships among all haplotypes identified from both ancient and modern deer of Hokkaido. Because the a-type is closest to all haplotypes from both ancient and modern deer, it is likely more ancestral among the haplotypes in Hokkaido.
Fig. 4 shows the neighbor-joining relationships among all haplotypes of Japanese sika deer (Nagata et al., 1998a; 1999) including new haplotypes identified in the present study. The three new haplotypes (g-, h- and i-types) were clustered together with six haplotypes (a- to f-types), and this cluster was supported with an 82% bootstrap value. The Hokkaido group was included in the northern group with the Honshu populations with a 100% bootstrap value. Meanwhile, sika deer of Yamaguchi, Kyushu, Tsushima, Yakushima and Kerama formed the southern group with a 99% bootstrap value.
Changes of geographical distribution patterns of mtDNA haplotypes between the past and present
The present study of ancient DNA indicates an evidence of changes of geographical distribution patterns in the Hokkaido sika deer population through bottleneck due to the heavy snow and hunting pressure in the beginning of Meiji Period. In winter of Hokkaido, plants (mainly bamboo grasses such as Sasa nipponia and Sasa senanensis) for sika deer's feeds are covered with the snow. To take these foods and escape much snow in winter, sika deer move into coniferous forests where the snow does not lie so much (Kaji, 1995). In Hokkaido, there are mainly three large coniferous forests: Daisetsu, Akan and Hidaka areas (Figs. 1a, 1b, 1c). Meanwhile, Nagata et al. (1998a; 1998b) found a low level of genetic variations in the modern sika deer population of Hokkaido. Based on the genetic data, they supposed that sika deer survived the heavy snow of 1879 and 1881 in these large forests and that the a-type increased in coniferous forests of the Akan area, b-type survived in forests of the Daisetsu area, and c-type survived in forests of the Hidaka area. Then, occupation of the three haplotypes could have expanded along the three forests, respectively.
By contrast, ancient DNA data of the present study demonstrated that distribution patterns of mtDNA haplo-types of archaeological deer were different from those of modern deer. Based on ancient DNA data and the report of Nagata et al. (1998a), we explain changes of population structures from the past to the present as follows. First, the a-type was major through huge areas of Hokkaido before the last bottleneck, that is, the heavy snow in 1879 and 1881. During those winters, a large number of the a-type survived and then expanded the distribution (Fig. 1a). Although distribution of the b-type in Ainu Culture Period was wider than the current one, this type was not major then. After the last bottleneck, the b-type survived around coniferous forests of the Daisetsu area, and its frequency in the local population became higher through genetic drift, and then expanded rapidly in northern parts of the Daisetsu area as well as coastal areas of Okhotsk Sea. In modern sika deer, the b-type is one of the most major types (Fig. 1b). Distribution of the c-type in Ainu Culture Period was wider than the current area and extended to eastern Hokkaido and Ishikari-lowland. However, the c-type disappeared in most of areas through the last bottleneck, and a small population with the c-type survived in coniferous forests of the Hidaka area (Fig. 1c). The other haplotypes such as d-, e-, and f-types were likely not major both before and after the last bottleneck.
Low genetic diversity in the Hokkaido population
Identification of the three new haplotypes (g-, h- and i-types) in the present study indicates that genetic diversity in the deer population before the bottleneck was greater than that of the modern deer population Hokkaido. The three new haplotypes were not major among ancient deer and identified restrictedly from three individuals at only Usu-oyakotsu Site (Table 3). Therefore, genetic diversity in Usu-oyakotsu Site of southern Hokkaido became the highest among archaeological sites. In modern sika deer, the e-type was found also in southern Hokkaido (Nagata et al., 1998a). These facts indicate that such minor haplotypes have been maintained in southern Hokkaido, where it is warmer and it snows less than in central and eastern Hokkaido.
The phylogenetic tree (Fig. 4) clearly supports that genetic diversity among the sika deer in Hokkaido has been much smaller than that in Honshu and Kyushu. The heavy snow in the beginning of Meiji Period is only a case of bottleneck in records (Kaji, 1995). Up to that time, there must have been many times of heavy snow which affected population structures of sika deer. In addition, environmental changes to affect vegetation could have occurred in Hokkaido. In such environmental changes, sika deer could have experienced overpopulation, crush or movement, resulting in changes of population structures. When subpopulation sizes were smaller, due to genetic drift, some haplotypes of mtDNA could have been extinct and others could have increased in number as mentioned above. In addition, genetic diversity of nuclear microsatellites in modern sika deer of Hokkaido is very low (Nagata et al., 1998b; Goodman et al., 2001). It is likely that three new types found in the present study disappeared through the bottleneck in Meiji Period.
Structure of tandem repeats in the mtDNA control region
The mtDNA control region of modern sika deer of Hokkiado commonly shared seven motifs of tandem repeats, each of which consist of approximately 40 bp (Nagata et al., 1998a, 1999) (Fig. 2). The unit number of tandem repeats differs among local populations on Japanese islands: four units at minimum (Yamaguchi, Tsushima, Goto islands, Yakushima) and seven at maximum (Hokkaido, Ciba, Siga, Gifu) (Nagata et al., 1999). The i-type having eight motifs identified from one archaeological deer has not been found from any modern deer, and it has two copies (duplicate) of the fifth unit, compared with the a-type. The sequence of this individual with the i-type was confirmed by multiple PCR amplifications and DNA sequencings. The finding of the i-type also supports that ancient deer were genetically more variable than the modern population.
The decline of genetic diversity was reported also by the nuclear microsatellite analysis of the modern sika deer in Hokkaido, compared with populations of Honshu (Nagata et al., 1998b). Therefore, nuclear DNA analysis is required to further investigate genetic diversity within the ancient population. However, because nuclear genes occur as only two copies per cell (vs. several thousand copies of mtDNA per cell), it is still difficult to PCR-amplify fragmented nuclear DNAs from archaeological bone samples. Actually we tried PCR amplification of nuclear DNAs in the present study, but it was impossible to obtain PCR products. Development of more efficient techniques is needed for ancient DNA analysis of nuclear genome of sika deer.
We thank S. Akaishi (Tomakomai City Museum), K. Fujita (Koshimizu Town Board of Education), A. Ishikawa (Kushiro Archaeological Research Center), Y. Masuda, I. Matsuda (Shire-toko Museum), K. Morioka, Y. Osada (Historical Museum of the Saru River) and N. Ohshima (Date City Board of Education) for supplying archaeological deer samples. We are also grateful to Prof. T. Nishimoto (National Museum of Japanese History) and H. Okada (Shiretoko Nature Foundation) for invaluable suggestions. This study was supported in part by Grants-in-Aid for Scientific Research and the 21st Century COE Program “Neo-Science of Natural History” at Hokkaido University from the Ministry of Education, Culture, Sports, Science and Technology, Japan, and a research grant from the Mitsubishi Foundation to R. M.
- H. Abe, N. Ishii, Y. Kaneko, K. Maeda, S. Miura, and M. Yoneda . 1994. A Pictorial Guide to the Mammals of Japan. Tokai Univ Press. Tokyo. in Japanese. Google Scholar
- Japanese Town Board of Education 1983. Hokkaido Akan Town Shimoninishibetsu Site. Akan Town Board of Education. Akan. in Japanese. Google Scholar
- Archaeological Exhibition Laboratory of Sapporo University 1998. Report of Excavation of Koshimizu Town Aoshimanai Site. Archaeological Exhibition Lab of Sapporo Univ. Sapporo. in Japanese. Google Scholar
- Biratori Town Board of Education 1997. Opaushinai 1 Site. Report of Cultural Assets of Biratori Town V. Biratori Town Board of Education. Biratori. in Japanese. Google Scholar
- Biratori Town Board of Education 1997. Penakori 1 Site. Report of Cultural Assets of Biratori Town VII. Biratori Town Board of Education. Biratori. in Japanese. Google Scholar
- Cooperation of Shiretoko Museum 1993. Report of Cultural Assets of Shari Town IV: Report of Excavation of Onnebetsu River Nishigawa-daichi Site. Cooperation of Shiretoko Museum. Shiretoko. in Japanese. Google Scholar
- Date City Board of Education 1993. Usu-oyakotsu Site amd Ponma Site. Date City Board of Education. Date. in Japanese. Google Scholar
- J. Felsenstein 1985. Confidence limits on phylogenies: An approach using the bootstrap. Evolution 39:783–791. Google Scholar
- S. J. Goodman, H. B. Tamate, R. Wilson, J. Nagata, S. Tatsuzawa, G. M. Swanson, J. M. Pemberton, and D. R. Mccullough . 2001. Bottlenecks, drift and differentiation: the population structure and demographic history of sika deer (Cervus nippon) in the Japanese archipelago. Mol Ecol 10:1357–1370. Google Scholar
- H. Harunari 2001. Relationship between the extinction of the big mammals and the human activities at the late Pleistocene in Japan. Bull Nat Mus Jpn His 90:1–52. in Japanese with English abstract. Google Scholar
- Hokkaido Government 1986. Result of a Survey Related to Sika Deer and Brown Bear Sightings in Hokkaido. Hokkaido Nature Preservation Division. Sapporo. in Japanese. Google Scholar
- Hokkaido Government 1994. Result of a Survey Related to Sika Deer and Brown Bear Sightings in Hokkaido. Hokkaido Institute of Environmental Sciences. Sapporo. in Japanese. Google Scholar
- T. Inukai 1952. On the deer in Hokkaido, its decline and coming back. Hoppou Bunka Kenkyu Houkoku (Studies from the Research Institute for Northern Culture of Hokkaido Univ) 7:1–45. in Japanese. Google Scholar
- K. Kaji 1995. Deer eruptions – A case study in Hokkaido, Japan. Honyurui Kagaku (Mammalian Science) 35:35–43. in Japanese. Google Scholar
- M. Kimura 1980. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide Ancient DNA of Sika Deer in Hokkaido 481 sequences. J Mol Evol 16:111–120. Google Scholar
- Koshimizu Town Office 1989. Furetoi Shell Mound. Koshimizu Town Office. Koshimizu. in Japanese. Google Scholar
- S. Kumar, K. Tamura, and M. Nei . 1993. MEGA: Molecular evolutionary genetic analysis. Version 1.01. Pennsylvania State Univ Park. Google Scholar
- Kushiro Archaeological Research Center 1996. Report of Kushiro City Nusamai Site III. Kushiro Archaeological Research Center. Kushiro. in Japanese. Google Scholar
- Kushiro Archaeological Research Center 1999. Report of Kushiro City Nusamai Site IV. Kushiro Archaeological Research Center. Kushiro. in Japanese. Google Scholar
- Kushiro City Museum 1975. Report of Kushiro City Katsurakoi Fushikokotan-chashi. Kushiro City Museum. Kushiro. in Japanese. Google Scholar
- R. Masuda, T. Amano, and H. Ono . 2001. Ancient DNA analysis of brown bear (Ursus arctos) remains from the archeological site of Rebun Island, Hokkaido, Japan. Zool Sci 18:741–751. Google Scholar
- J. Nagata, R. Masuda, K. Kaji, M. Kaneko, and M. C. Yoshida . 1998a. Genetic variation and population structure of the Japanese sika deer (Cervus Nippon) in Hokkaido Island, based on mitochondrial D-loop sequences. Mol Ecol 7:871–877. Google Scholar
- J. Nagata, R. Masuda, K. Kaji, K. Ochiai, M. Asada, and M. C. Yoshida . 1998b. Microsatellite DNA variations of sika deer, Cervus Nippon, in Hokkaido and Chiba. Mammal Study 23:95–101. Google Scholar
- J. Nagata, R. Masuda, B. H. Tamate, S. Hamasaki, K. Ochiai, M. Asada, S. Tatsuzawa, K. Suda, H. Tado, and M. C. Yoshida . 1999. Two genetically distinct lineages of the sika deer, Cervus nippon, in Japanese islands: Comparison of mitochondrial D-loop rejoin sequences. Mol Phylogenet Evol 13:511–519. Google Scholar
- N. Saitou and M. Nei . 1987. The neighbor-joining method: A new method for reconstructing phylogenetic trees. Mol Biol Evol 4:406–425. Google Scholar
- J. Sambrook, E. F. Fritsch, and T. Maniatis . 1989. Molecular Cloning: A Laboratory Manual. 2nd ed. Cold Spring Harbor Lab Press. NY. Google Scholar
- Tomakomai City Archaeological Research Center 1987. Benten Shell Mound I -Report of Excavation of Ainu Shell Mound after the End of Edo Period-. Tomakomai City Archaeological Research Center. Tomakomai. in Japanese. Google Scholar
- Tomakomai City Archaeological Research Center 1988. Benten Shell Mound II -Report of Excavation of Ainu Shell Mound after the End of Edo Period-. Tomakomai City Archaeological Research Center. Tomakomai. in Japanese. Google Scholar
- Tomakomai City Archaeological Research Center 1989. Benten Shell Mound III -Report of Excavation of Ainu Shell Mound after the End of Edo Period-. Tomakomai City Archaeological Research Center. Tomakomai. in Japanese. Google Scholar
- Tomakomai City Board of Education 1975. Report of Excavation of “Deer Meat Factory in Reclaimed Bibi” in Bibi Tomakomai City -Report of Cultural Assets of Tomakomai City-I. Tomakomai City Board of Education. Tomakomai. in Japanese. Google Scholar
- T. Watanobe, N. Ishiguro, N. Okumura, M. Nakano, A. Matsui, H. Hongo, and H. Ushiro . 2001. Ancient mitochondrial DNA reveals the origin of Sus scrofa from Rebun Island, Japan. J Mol Evol 52:281–289. Google Scholar