Genetic population structure of invasive raccoons (Procyon lotor) in Hokkaido, Japan: Unique phenomenon caused by pet escape or abandonment.

Phylogeographic studies can resolve relationships between genetic population structure of organisms and geographical distributions. Raccoons have become feral in Japan, and in Hokkaido island, they have been rapidly increasing in number and spreading since the 1970s. We analyzed mitochondrial (mtDNA) and microsatellite DNA to understand the current phylogenetic distribution and invasive founder events. Overall, Hokkaido raccoons maintained high genetic diversity (i.e., the level of heterozygosity was comparable to the original habitat, North America). Based on mtDNA distribution and microsatellite diversity, Hokkaido raccoons were divided into six management units. However, mtDNA haplotype distributions and genetic structures based on microsatellites did not always correspond to each other (e.g., two geographically and genetically separated populations showed similar mtDNA distributions). In addition, a high degree of genetic admixture was observed in every unit, and the degree of genetic differentiation was low even between regions separated by long distances. Compared with other countries in Europe where genetic distribution of introduced raccoons is more clearly structured, the current results represent a unique and complex phenomenon of pet escape/abandonment in Hokkaido: i.e., genetically related colonies were introduced into multiple regions as founder events, resulting in the current state in which raccoons are not clearly genetically differentiated even 40 years after introduction.


Results
Mitochondrial diversity of raccoons in Hokkaido and other regions. With 392-bp fragment, seven haplotypes were detected with 17 polymorphic sites (Table 1): NA (n = 234), HT (n = 117), KK (n = 92), AS (n = 14), AB (n = 62), AH (n = 5), and IW (n = 2) in 526 raccoons collected from 44 municipalities in ten subprefectures in Hokkaido. In a longer 709-bp fragment, an additional nine polymorphic sites were detected, and every sample from the same haplotype showed an identical sequence. Nucleotide sequence data are available in the DNA Data Bank of Japan database as well as GenBank with the following accession numbers: IW: LC455749, HT: LC455750, KK: LC455748, AS: LC455753, NA: LC455747, AB: LC455752, and AH: LC455751. Some haplotypes detected in this study were identified in reports from Japan 26,27 and North America 18 (Table 1)
fragment length is different among the studies (Fig. 1). Polymorphic sites among haplotypes IW, HT, and KK were detected only in base positions 58 and 160, which were not examined in Cullingham's study 18 . Therefore, these three haplotypes were shown to be the same as PLO32 in North America (Table 1). Haplotype NA was found in almost all subprefectures in the sampling area (Fig. 2B) and detected in about half of the samples (44.5%) investigated in this study. All of the individuals in northern Hokkaido were identified as this haplotype. Haplotype HT was found in about half of the municipalities but mainly in southern-central Hokkaido (Fig. 2B). The majority (95%) of individuals in Tokachi subprefecture were classified as this haplotype. Haplotype KK was found mainly in west-central Hokkaido (Fig. 2B). Haplotypes AS, IW, and AB were found only in limited municipalities (Fig. 2B). Haplotype AH was found only in five individuals in a municipality of Kamikawa Subprefecture (Fig. 2B). This haplotype differed by one nucleotide in base position 227 compared with haplotype AB. The largest number of haplotypes were detected in Kamikawa and the east Sorachi area (six haplotypes), with the second largest in Ishikari and the southwest Sorachi area (five haplotypes).
Locus-specific diversity in microsatellite analysis. A total of 326 raccoons, which were randomly chosen with maximum 10 samples from each 44 municipality, were genotyped at ten microsatellite loci. The number of alleles per locus ranged from 6 to 19. Single locus heterozygosities observed (H o ) ranged from 0.66 to 0.87, with 0.75 on average (0.07 standard deviation (SD)). The expected heterozygosities (H e ) ranged from 0.65 to 0.89, with 0.78 on average (0.08 SD) ( Table 2). HWE tests and LD tests were conducted using the Markov chain method (1000 dememorization steps, 1000 batches, and 1000 iterations per batch). Significant values in HWE were different for one locus from five subprefectures (P < 0.05) but were from only two subprefectures following a strict Bonferroni correction for multiple tests (P = 0.05/10 = 0.005); only one locus deviated from HWE in the two subpopulations. Out of 45 pairs of loci in LD tests, one to three pairs were significantly linked by Bonferroni correction in three subprefectures (P = 0.05/45 = 0.0011). However, there were no loci among any pairs that were consistently significant over past studies 18,[28][29][30] , and no value was significant beyond several populations in this study. There was no consistent evidence for LD or null alleles across any populations. Therefore, all loci were included in further analyses. population structure and genetic diversity. Structure Harvester results analyzed based on 44 municipalities showed the lowest ΔK when K = 4 (Supplement file 1). However, log-likelihoods were not high at K = 4 and increased beyond K = 4. The highest values were obtained for K = 5, and, therefore, we chose K = 5 as an appropriate number of clusters. In the STRUCURE analysis, 44 municipalities were divided into five clusters, which consisted of the following subprefectures: 1: Soya and Rumoi; 2: Kamikawa, Sorachi, and Ishikari; 3: Shiribeshi and Iburi; 4: Hidaka; 5: Tokachi and Doto (Fig. 2C). However, mtDNA distributions in the area of Kamikawa and east Sorachi were distinct from those in the area of southwest Sorachi and Ishikari. Therefore, based on the mtDNA distributions and microsatellite diversity with K = 5, we divided raccoon groups into six management units in Hokkaido (Fig. 2D), as follows: Unit1 (Soya and Rumoi subprefectures), Unit2 (Kamikawa and east Sorachi subprefectures), Unit3 (the region covering southwest Sorachi, Ishikari, and a part of Iburi subprefectures), Unit4 (Shiribeshi, west Iburi, and west Ishikari subprefectures), Unit5 (Hidaka and east Iburi subprefectures), and Unit6 (Tokachi Subprefecture). With 321 raccoons forming these six management units, H o ranged from 0.726 to 0.784 with an average of 0.75 (0.02 SD). H e ranged from 0.731 to 0.787 with an average of 0.761 (0.02 SD) ( Table 3). Pairwise F ST values among the six units were statistically differentiated between any two units, and values were low, ranging between 0.0045 to 0.0667 (Table 4).

Discussion
In the present study, mitochondrial and nuclear diversity was examined in feral raccoons introduced into Hokkaido, Japan. This is the first report to study raccoon population dynamics based on genetic differentiation in Japan and to compare genetic characteristics among this animal in Japan with an original habitat in North America.
The 662-bp sequences detected in this study covered fragments identified in previous studies, including 289 bp in Japan 27 , 268 bp in North America and West Indian islands 16,17 , and also 467 bp in eastern North America 18 (Fig. 1). The sequences also partially covered 521 bp out of a 550 bp fragment from Germany 23 . Here, seven haplotypes were detected from 526 samples. Five of these were reported previously in central Hokkaido 27 . When adding RMT06, which had been found only in the study by Takada-Matsuzaki et al. 27 but not in the present study, founder events in Hokkaido were established from at least eight female haplotypes ( Fig. 2A). This could be an underestimation, given that founder phylogenetic diversity can diminish during changing demographic periods based on the number of founder females and survivorship. Compared with reports in native habitats in Northeastern America, where 76 haplotypes were detected from 311 samples 18 , the number of haplotypes in Hokkaido is small. Such low mitochondrial diversity was also reported in other introduced areas in Europe, with two haplotypes identified from 58 samples in Spain 13 , six haplotypes from 193 samples in Germany 23 , and four haplotypes from 72 samples from the border area in Poland, Czech Republic, and Germany 22 . Rapid expansion from small populations can involve serial bottlenecking with progressive loss of allelic diversity, so that populations in more recently colonized locations could contain less genetic diversity 31 .
Four haplotypes out of eight identified in this study were also detected in North America 18 and five haplotypes out of the eight were also detected in the other locations within Japan 32 , indicating that raccoons from the same origin group in North America were introduced at several remote places in Japan as founder events. Another possibility is that re-introduction (as a second release in the prefectures) from the initial introduced regions into the non-habitat area, happened after the first founder event as discussed by Takada-Matsuzaki et al. 27 Haplotype NA matched haplotype PLO2, which was mentioned in Cullingham's study 18 . PLO2 is the most common raccoon haplotype, is widespread in North America, and has also been reported in non-native areas, such as Germany, Spain, Poland, and Czech Republic 13,22,23 . Haplotype NA was detected in introduced areas across countries, which strongly suggests that the same genetic population of this haplotype was managed under captivity from wild populations in North America, exported as pets or as a fur-bearing animal and was introduced into non-native regions.
In this study, 326 individuals were divided into six groups based on mtDNA polymorphism, STRUCTURE analysis, and habitat regions, by use of the municipality as a unit. We note that this classification might not be biologically relevant for wild animals that move regardless of the frontiers between municipalities. In Hokkaido, however, most eradication programs of raccoons have been implemented by each municipality. Therefore, this Hokkaido (n = 326) Chicago (n = 323) 28 Table 2. Locus-specific diversity measures at ten microsatellite loci for 326 raccoons in Hokkaido, Japan, and other reports from native raccoon habitat in North America. The number of alleles (Na), observed heterozygosity (Ho), expected heterozygosity (He). The data set from Chicago was reported in Santonastaso et al. 28 . The data set from Illinois was reported in Hauver et al. 33 .   study focused more on classifying management units, rather than on finding the ecological boundary between sub-populations. Genetic diversity at these ten microsatellite loci in Hokkaido raccoons is low in Na though it was at the same level in H o and H e when compared to that in a native habitat of North America (in Chicago 28,33 and Illinois 34 ) ( Table 2). When genetic variability defined by clusters was evaluated in the six management units (Table 3), Ar, H o, and H e were high compared to Germany where feral raccoons are widespread throughout the country 24 . Values were low in six clusters of 407 feral German raccoons: Ar: 3.6-4.9, H o : 0.52-0.64, and H e : 0.53-0.65. These differences demonstrate that nuclear genetic diversity in Hokkaido has stayed polymorphic at the same level as native habitat, in contrast to other invaded habitats.

Number of samples
In other raccoon-invaded areas, genetic diversity at the nuclear level in Spain was lower than in Missouri, a native area 13 . Raccoons were first reported in Spain relatively recently, in 2003. The microsatellite analysis consisted of 58 samples and showed only two mtDNA haplotypes captured in the limited area of Madrid. Therefore, diversity may be lower than observed in our study and in reports from the United States. In contrast, raccoons have been in Germany much longer than in Spain. There, raccoons were introduced in the 1930s, but the populations were derived from two separate founding events and one additional escape event 24 . These small numbers of founding events likely led to low polymorphisms seen in Germany 24 . The high genetic polymorphisms observed in this study indicate that the founder event in Hokkaido was relatively complex compared to that in Spain or Germany. The mtDNA of raccoons in Hokkaido showed allopatric, whereas the genome DNA showed very low pairwise F ST values between any two management units regardless of the distance of the two units. Almost all the F ST values (0.005-0.067) in the current study was lower than the lowest value observed in Germany (0.049) 23 . This indicates that several maternal populations, with different genetic backgrounds, were introduced into different regions in Hokkaido. However, these original populations were closely related to each other, resulting in limited differentiation among Units even 40 years after raccoons were introduced.
In Japan, raccoons were mainly imported as pets and introduced into the wild afterwards 35 . To our best knowledge, there are no documents recording when and how many raccoons were imported into Japan 7 , though the results of this study reflect that founder populations in Hokkaido likely derived from genetically related individuals from a captive breeder in North America or the main island of Japan. When allelic richness was compared among the six management units in Hokkaido, Unit2 and Unit3 showed relatively high genetic diversity and also highly diverse mtDNA (Table 3). Unit3, covering Sorachi, Ishikari, and Iburi subprefectures (the most human populated regions), which includes Sapporo, the largest city in Hokkaido. Unit2 has the second largest city, Asahikawa, suggesting that raccoons from multiple haplotypes were released into highly populated regions. Historical evidence demonstrates that Ishikari, Sorachi, and Kamikawa subprefectures were regions where feral raccoons or their traces were detected in the 1990s, earlier than in other regions (Hokkaido governmental reports and Ikeda 7 ). A report showing that national parks with high surrounding human population density had significantly more alien species 36 also supports our findings. Additionally, in Hokkaido, genetic polymorphism was varied but conserved, not only in highly populated areas, but also in regions throughout the island. This result suggests that several founder females with different haplotypes were introduced sympatrically and/or allopatrically, and cross-municipality migrating has resulted in the genetic diversity observed in Hokkaido.
Ikeda 6 reported that the naturalization of raccoons in Hokkaido was due to multiple escapes and abandonment of pets. They also reported that additional human-induced transport and secondary releases have allowed habitat expansion in more remote locations within Hokkaido. Raccoons tend to disperse along the rivers and water passages 37 . Females are often philopatric, while males are more likely to disperse 37,38 . In the native habitat, the majority of raccoons relocate less than 5 km 39 , meanwhile they sometimes disperse more than 20 km in the absence of major landscape barriers 30 . Genetic analysis here suggests that geographical barriers prevent migration of feral individuals and reduce genetic admixture between some populations. The mtDNA haplotype of Unit1 and a part of Unit4 were the same; genetic structure in microsatellite clusters was different between these two units, indicating that there was less genetic intercrossing between north and southwest regions in Hokkaido. Geographical distances between these regions also support that there has been less genetic admixture between these units after individuals were introduced from the same origin. Unit2 was limited to Kamikawa and Sorachi Subprefectures, and Unit6 in Tokachi Subprefecture, which are neighboring regions. Although mtDNA haplotypes were different between these two units, STRUCTURE analysis inferred possible genetic admixture, suggesting that male individuals might have been dispersed among these regions.
Unit5 in Hidaka Subprefecture and Unit6 are also neighboring regions and have the same mtDNA haplotype constitution. However, STRUCTURE results indicated that there was less individual contact. There is approximately 150 km of the Hidaka mountain range between Hidaka and Tokachi, which could have prevented genetic mixture by serving as a geographical barrier, resulting in each population having maintained their own genetic structure. It is thus important for neighboring areas to cooperate with each other to control the spread of invasive raccoons. In the case of Hidaka and Tokachi, there could be effective management considering the geographical barriers 19 . Although human-induced re-introduction and transport were reported at the beginning of raccoon naturalization, the Invasive Alien Species Act of 2005 strictly regulates transportation of living animals. In this situation, biological information, such as geographical isolation or populational genetics, will be helpful for effective eradication and control 19 .

Summary
1. In Hokkaido, founder events of raccoon introduction occurred sympatrically and allopatrically in separate regions, which reflects a unique phenomenon of pet escape/abandonment with different maternal populations (consisting of at least eight mtDNA haplotypes) that were closely related with each other. (2020) 10:8108 | https://doi.org/10.1038/s41598-020-64526-y www.nature.com/scientificreports www.nature.com/scientificreports/ 2. Genetic diversity of raccoon has stayed polymorphic at almost same level as native habitat, and a high degree of genetic admixture was detected in regions where humans were highly populated. 3. Some populations have maintained their own genetic structure, in areas of geographical barriers.
Compared with other countries in Europe where genetic distribution is more clearly structured in introduced raccoons, the current results show that genetically related colonies were introduced into multiple regions as founder events, resulting in the current state where raccoons are not clearly genetically differentiated even 40 years after introduction.

Materials and Methods
Sample collection and DNA extraction. Hair To confirm amplification, 5 µL of the products were electrophoresed on a 1.5% agarose gel. The products were purified using a FastGene Gel/PCR Extraction Kit (NIPPON Genetics Co., Ltd.). The products were sequenced using a Big Dye Terminator v1.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA). The sequences were analyzed using an ABI PRISM 310 Genetic Analyzer (Applied Biosystems) and differentiated into several haplotypes. Subsequently, to compare these 392-bp sequences in Hokkaido with those in North America 18 , another reverse primer, ProcR-add (5′-CCATGAATTAAACTGCACCA-3′), was designed for amplification of longer 709-bp fragments based on the sequence reported by Frantz et al. 23 . In 46 samples, which were picked from seven haplotype groups, covering all detected subprefectures, longer 709-bp fragments were amplified with primers Proc L and ProcR-add and sequenced as described above, in order to confirm that each sequence detected in Hokkaido was identical to the haplotype reported in North America (Fig. 1).

Molecular data analysis.
Seven haplotypes of the D-loop fragment sequence and one haplotype reported in a previous study 26 were aligned and a neighbor-joining phylogenic tree was constructed using MEGA5 software (http://www.megasoftware.net/) 40 . Genetic distances, based on alignment gaps among the haplotypes, were estimated from all substitutions and by using the two-parameter distance method 41 .
Microsatellite genotyping. A total of 326 individuals from all 44 municipalities were genotyped at 10 microsatellite loci (Table 5). From municipalities with>10 samples, 10 samples were randomly selected. From the other municipalities, with 10 or less samples, all samples were used for the microsatellite analysis. Primers and primer combinations used in multiplex PCR assays are shown in Table 5. Multiplex PCR reactions were performed in a total volume of 15 µL, consisting of 1 µL of DNA solution, 0.075 µL of Kit Mix 1 (Multiplex Assay Kit; Takara Bio Inc., Shiga, Japan), 7.5 µL of Kit Mix 2, 0.5 µL of primer mix (0.25 µM each), and 5.925 µL of PCR-grade water. The mixture was heated to 94 °C for 1 min, followed by 40 cycles of 30 s at 94 °C, 1 min at 52-56 °C, and 1 min at 72 °C, with final 10 min at 72 °C. One microliter of product was diluted in 120 µL of distilled water, and 1 µL of this diluted solution was mixed with 0.125 µL of GeneScan 500 LIZ size standard (Life Technologies Japan Ltd., Tokyo, Japan) and 10 µL of Hi-Di formamide (Life Technologies Japan Ltd.). After denaturation at 95 °C for 3 min and cooling on ice, each sample was analyzed with the ABI PRISM 310 Genetic Analyzer (Life Technologies Japan Ltd.). Allele size was determined using GeneScan ver. 4.1 software (Life Technologies Japan Ltd.). When PCR amplification was weak for any locus, a single PCR reaction was performed with the same primers. Data analysis. The software CERVUS 3.0.7 42 was used to calculate loci characteristics, including number of alleles (Na), observed and expected heterozygosities (H o and H e ), and polymorphic information content. In addition, we tested whether there was any evidence suggesting a departure from Hardy-Weinberg equilibrium (HWE) of ten microsatellite loci and linkage disequilibrium (LD) between populations using the web-based program GENEPOP 4.2 43,44 .