Post-bottleneck increase in mitochondrial DNA diversity in Yaku sika deer (Cervus nippon yakushimae) on Yakushima Island, Japan

TL;DR

This study examines how a population bottleneck followed by genetic diversification shaped the mitochondrial DNA of Yaku sika deer on Yakushima Island.

Contribution

The study reveals a post-bottleneck increase in genetic diversity and historical dispersal events in Yaku sika deer.

Findings

Yakushima's Yaku sika deer show 18 haplotypes with regional variation and two subpopulation structures.

A historical dispersal event is inferred to have occurred approximately 4,900–6,500 years ago after a population bottleneck.

The Yakushima and Tanegashima populations are genetically distinct from other southwestern Japanese populations.

Abstract

We investigated the genetic diversity and divergence of sika deer (Cervus nippon) descended from survivors of a catastrophic pyroclastic flow event approximately 7,300 years ago by analyzing the 894 base pairs mitochondrial DNA control region extracted from fecal samples. We focused on Yaku sika deer (C. n. yakushimae) across the remote islands Yakushima and Kuchinoerabujima in southern Japan. On Yakushima, 18 haplotypes were detected, with significant regional variation. Two subpopulation structures on Yakushima aligned with the distribution of pyroclastic flow deposits. Two haplotypes were found on Kuchinoerabujima, one shared with Yakushima. Reconstructed haplotype networks with sika deer from various regions of the Japanese Archipelago suggested that a historical simultaneous dispersal event occurred in the Yakushima population. The Yakushima and Tanegashima populations (C. n. mageshimae) formed their own clusters and were clearly separated from other southwestern Japanese populations. The lack of common haplotypes between these populations suggests that the Yakushima population experienced sufficient generational changes for genetic drift to fix region-specific mutations after geographical isolation approximately 37,500 years ago. The high genetic diversity of the current Yakushima population may have been acquired through simultaneous dispersal approximately 4,900–6,500 years ago after the pyroclastic flow-induced population bottleneck approximately 7,300 years ago.

The online version contains supplementary material available at 10.1038/s41598-025-23191-9.

Figures

Study site, sample collection points, and area divisions used in our analysis. (a) Map of the Japanese Archipelago and the locality used for haplotype network analysis. The star (★) represents the Kikai Caldera. An enlarged view of the area enclosed by the solid line is shown in “b”. (b) Enlarged map of Yakushima and Kuchinoerabujima. Black circles (●) represent the sample collection points. The black triangles (▲) represent the highest elevation points on each island. The A–I areas used in the analysis are circled by solid lines. This Figure was created based on the GSI Map Vector of the Geos

Haplotype network based on the mitochondrial DNA control region of sika deer (*Cervus nippon*) in the Japanese Archipelago. We analyzed 110 haplotypes (Supplementary Table S1) by combing data from the NCBI registry and this study. Yaku sika deer haplotypes are enclosed in dashed boxes and enlarged Figures are shown in a separate box “a”; Yk01–Yk18 and Kc01 represent the haplotype names in this study, and W represents the sequence ([AB218689](AB218689)) from Wada et al.^[40](#CR40)^. The dashes in the branches represent mutation points and the small black circles represent hypothetical haplotyp

Genetic population structure of Yaku sika deer (*Cervus nippon yakushimae*) and thickness of the pyroclastic flow deposited on Yakushima. Black triangle (▲) represents the highest elevation point (1,936 m asl). Thickness of pyroclastic flow deposits estimated by spatial interpolation (Gaussian model, 5 km mesh) based on Geshi^[36](#CR36)^. SAMOVA ver. 2.0^[92](#CR92)^ was used to divide the Yaku sika deer into two genetic populations, i.e., north (open circles ○; areas B/C/G/I) and the other areas (filled circles ●; areas A/D/E/F/H: Fig. [1](#Fig1)b). However, Area A, surrounded by the dashed

Mismatch distribution of Yaku sika deer (*Cervus nippon yakushimae*) based on the mitochondrial DNA control region 894 base pairs. The thick solid line represents the observed mismatch distribution (—●—). The thin solid line represents the demographic expansion model (—□— ). Two dashed lines indicate 95% confidence intervals (----).

Maximum likelihood tree based on the mitochondrial DNA control region for the 110 haplotypes of the sika deer (*Cervus nippon*) in Japan with Wapiti (*C. canadensis*) as the outgroup. Numbers of nucleotide substitutions per site indicated by scale are the Tamura 3-parameter distances^[97](#CR97)^. Numbers on branches indicate bootstrap values after 1,000 replicates. See Supplementary Table S4 for haplotype information.

Full text

Introduction

Genetic diversity is one of the most important factors for the conservation of wild animal populations because low genetic diversity leads to decreased population fitness^1^ and may even result in vulnerability to extinction by factors such as inbreeding and loss of adaptability to environmental changes^2,3^. Population bottlenecks, i.e., extreme reductions in population size, have occurred in various populations, often leading to a decrease in genetic diversity^4–6^. Understanding the recovery process of genetic diversity after bottlenecks will help conserve populations and provide insights on the formation of local genetic structures and species speciation. In this study, we analyzed the genetic diversity of an isolated population of sika deer (Cervus nippon) that likely experienced a bottleneck due to catastrophic events and examined when and how its genetic diversity was acquired.

Sika deer are widely distributed in East Asia, extending from central China (western end) to Japan (eastern end) and from the easternmost part of Russia (the northern end) to southern China and Vietnam (southern end); however, most of the continent’s wild populations are extinct or threatened today^7^. In contrast, sika deer populations in the Japanese Archipelago temporarily declined due to habitat destruction during the modern and contemporary periods, but have been recovering since the 1970s^8–10^. The sika deer populations across the Japanese Archipelago are classified into six subspecies in each habitat (Fig. 1a,b): C. n. yesoensis (Hokkaido; 77,984 km^2^), C. n. centralis (Honshu; 227,943 km^2^ and Tsushima Islands), C. n. nippon (Shikoku; 18,298 km^2^, Mainland Kyushu; 36,783 km^2^ and Goto Islands), C. n. mageshimae (Tanegashima; 445 km^2^ and Mageshima; 8 km^2^), C. n. yakushimae (Yakushima; 504 km^2^ and Kuchinoerabujima; 36 km^2^), and C. n. keramae (Kerama Islands)^11^. However, studies based on mitochondrial DNA (mtDNA)^12–14^ and nuclear DNA microsatellite analyses^15^ have reported the existence of two distinct lineages—northern and southern—highlighting the need to elucidate genetic structure as a basis for defining conservation management units^16–18^. These lineage boundaries are consistently recognized in Southern Honshu and Shikoku^12–15^, with divergence times estimated at approximately 0.35 million years ago (mya) based on mtDNA D-loop partial sequences^12^, 0.41 mya based on complete cytochrome b sequences^19^, and 0.53 mya based on mitogenomic data^20^.Furthermore, sika deer populations in Japan exhibit karyotype polymorphisms (2n = 64, 66, or 68), and their hybrids (2n = 65 and 67) have been identified^21–23^. According to Harada^21^, the northern lineage has 2n = 68 chromosomes, which is the oldest and most widely distributed on Hokkaido, Northern Honshu, most of Central Honshu, and Shikoku (Fig. 1a). For the southern lineage, most Southern Honshu populations have 2n = 66 chromosomes, while most Kyushu populations have 2n = 64 chromosomes.

Fig. 1. Study site, sample collection points, and area divisions used in our analysis. (a) Map of the Japanese Archipelago and the locality used for haplotype network analysis. The star (★) represents the Kikai Caldera. An enlarged view of the area enclosed by the solid line is shown in “b”. (b) Enlarged map of Yakushima and Kuchinoerabujima. Black circles (●) represent the sample collection points. The black triangles (▲) represent the highest elevation points on each island. The A–I areas used in the analysis are circled by solid lines. This Figure was created based on the GSI Map Vector of the Geospatial Information Authority of Japan (GSI) (https://maps.gsi.go.jp/vector/#7/36.104611/140.084556/&ls=vstd&disp=1), using the GSI Map Sheet Tool (Ver. 2.0.0.1) provided by GSI.

The Yaku sika deer (C. n. yakushimae) examined in this study is a subspecies of sika deer and is distributed exclusively on Yakushima and Kuchinoerabujima Islands^24^. Although the Yaku sika population was clustered with the southern Japanese lineage based on mtDNA^14,25^ and nuclear DNA phylogeny^15^, the genetic distance was significantly different from those of the other populations of the southern Japanese lineage^14,15,25,26^. Polymorphisms of Y chromosome simple sequence repeats (SSRs) was also clearly differentiated from those of the other Japanese populations^18,20^. The Yakushima populations have 2n = 68 chromosomes, which is similar to that of the northern lineage populations. Thus, the karyotype characteristics of the Yakushima population could not be classified simply based on geographical location. The Yaku sika deer is the smallest among the subspecies of the Japanese sika deer and has shorter limbs; these morphological traits are thought to be due to genetic characteristics^20,26–28^. These deer are thought to have evolved in Yakushima, an island that is not inhabited by medium- and large-sized carnivorous predators^29^. Miniaturization and limb shortening are typical morphological traits that frequently evolved in ungulates inhabiting small islands without predators^30^. In addition to these morphological features, the ecology^31–33^ and life history^34^ of these deer clearly differed from those of other sika deer subspecies. In the western lowland of Yakushima, females give birth for the first time at 3–5 years of age, and three-pointed antlers appear in males aged 4–5 years or older (Agetsuma and Agetsuma-Yanagihara unpublished data). The average body weight of females and males older than 4 years is 21 kg and 28 kg, respectively (Agetsuma and Agetsuma-Yanagihara unpublished data). The mean annual home range size (expressed as a 90% fixed kernel) for females and males aged more than 4 years is 12 ha and 36 ha, respectively, which is smaller than that of other subspecies^35^.

Approximately 7,300 years ago, the Yaku sika deer population would have encountered catastrophic damage due to the eruption of a submarine volcano named the Kikai Caldera (Fig. 1a,b). Pyroclastic flow from the Kikai Caldera (Koya pyroclastic flow) ran across the islands and deposited sediments over most parts of Yakushima and Kuchinoerabujima^36^. Wildlife populations in these islands were greatly impacted by the pyroclastic flow. Although migratory individuals from different populations generally bring diversity to genetically impoverished populations, this requires connectivity between habitats that allows for movement between populations. However, for populations on isolated islands, recovery of genetic diversity is more difficult. In fact, the sympatric Yakushima macaques (Macaca fuscata yakui) is thought to have experienced a strong bottleneck after the catastrophe^37^. Thus, sympatric macaques and Yaku sika deer might have been similarly impacted. Nonetheless, Terada et al.^17^ performed a microsatellite genetic structure analysis and showed that Yaku sika deer have high genetic diversity despite the small island population and are separated into two or four subpopulations (northern and southern, or northern, southern, eastern, and western parts of Yakushima); however, sufficient samples from the wildlife protection areas and mountainous areas were not analyzed. We examined whether the Yaku sika deer experienced a bottleneck, and if so, how and when the current genetic diversity was acquired. Our analyses will help reveal how genetic structure and diversity are reformed in isolated populations. mtDNA is a better marker for detecting genetic variation within a single species than nuclear DNA because it does not involve recombination and mutations accumulate rapidly due to the high rate of base substitution^38,39^. Thus, mtDNA is suitable for investigating regional differences within small islands such as Yakushima. Therefore, by examining the mtDNA of deer from whole areas of Yakushima Island, we investigated the process by which genetic diversity was obtained in the island population.

Results

Genetic relationships among populations from different areas of Yakushima

We determined the partial mtDNA control region sequence in 380 out of 419 samples in Yakushima and 2 out of 5 samples in Kuchinoerabujima. A total of 303 samples from Yakushima and 1 sample from Kuchinoerabujima yielded more than 99% of the 894 base pairs (bp) partial sequence, beginning at the 48th bp of the reference sequence registered in the NCBI database (Accession No. AB279718.1; https://www.ncbi.nlm.nih.gov). This 894 bp partial sequence of the mtDNA control region in Yaku sika deer contained 13 single-nucleotide substitutions (all transitions) and two single-nucleotide length polymorphisms. By combining these data, 19 haplotypes were identified from Yakushima and Kuchinoerabujima (Supplementary Table S1 online). The other sample from Kuchinoerabujima, Kc02, yielded a slightly shorter sequence of 841 bp, but it was identical to the Yakushima haplotype Yk06 within the aligned region (Supplementary Table S1 online). A total of 20 haplotypes, including Kc02, were deposited in the DNA Data Bank of Japan (DDBJ; https://www.ddbj.nig.ac.jp/index.html) under Accession Nos. LC664029–LC664048. Due to its shorter sequence length, Kc02 was excluded from subsequent analyses. A haplotype network constructed using the haplotypes obtained in this study, together with the sequence registered in NCBI (AB218689) from Wada et al.^40^, exhibited a radial shape centered on Yk01, with most haplotypes connected by a single mutation (Fig. 2a).

Fig. 2. Haplotype network based on the mitochondrial DNA control region of sika deer (Cervus nippon) in the Japanese Archipelago. We analyzed 110 haplotypes (Supplementary Table S1) by combing data from the NCBI registry and this study. Yaku sika deer haplotypes are enclosed in dashed boxes and enlarged Figures are shown in a separate box “a”; Yk01–Yk18 and Kc01 represent the haplotype names in this study, and W represents the sequence (AB218689) from Wada et al.^40^. The dashes in the branches represent mutation points and the small black circles represent hypothetical haplotypes. Created using Network 10.2.0.0 (https://www.fluxus-engineering.com/sharenet.htm).

Among the 303 samples from Yakushima, the total haplotype diversity (h) was 0.9006, with the diversity from each area in Yakushima ranging from 0.6931 (Area H) to 0.8480 (Area F). The total nucleotide diversity (π) was 0.00293, with the diversity from each area ranging from 0.00063 (Area H) to 0.00330 (Area C) (Table 1). Yk01, Yk05, Yk06, Yk07, and Yk11 were distributed over a wide area of the island, Yk09 was distributed only in the central to northern part of the island (areas B, C, G, and I), while Yk04 was primarily (96.5%) distributed in Area A (Supplementary Table S2 online).

Table 1 Summary statistics of nucleotide sequence data calculated based on the mitochondrial DNA control region 894 base pairs of Yaku sika deer (Cervus nippon yakushimae) in Yakushima calculated using Arlequin ver 3.5.2^93^.AreaNo. of samplesNNo. of polymorphic sitesSNo. of haplotypesHHaplotype diversityhNucleotide diversityπMean No. of differenceskA1041080.69530.002772.3373B18640.73860.002662.2422C22970.83550.003302.7797D30560.82990.001941.6349E21660.76670.001251.0526F19670.84800.001671.4090G38860.79520.002402.0279H28560.69310.000630.5313I23660.69960.002181.8337Northern: B/C/G/I20213150.86440.002472.0856Southern: A/D/E/F/H10112110.83230.002662.2465Total30314180.90060.002932.4732

With respect to the relationship between geographical location and genetic differentiation, spatial analysis of molecular variance (SAMOVA) revealed a fixed index FCT = 0.2250 (p = 0.00587) when the number of groups K equaled 2 (i.e., northern [Areas B, C, G, and I] and southern populations [Areas A, D, E, F, and H]). At K = 3, Area A was distinguished from the southern population and became a single group. At K = 4 and 5, Area A remained the lone group. The northern population was located on the thick sediment of the Koya pyroclastic flow, whereas the southern population (excluding Area A) was on the thin sediment (Fig. 3). Among areas of the southern population, Area A had exceptionally thick sediment. Genetic differentiation between northern and southern populations examined by analyses of molecular variance (AMOVA) indicated that 26.27% of the variation in the total sample population was due to the differences between areas (FST = 0.26267, p < 0.0001, Supplementary Table S3 online). In addition, pairwise FST values between areas showed significant differentiation among many areas. Especially, the FST of the population in Area A was significantly different from that of all the other eight areas (p < 0.001; Table 2).

Table 2. Degree of genetic differentiation among nine areas (A–I) by pairwise FST values based on the mitochondrial DNA control region 894 base pairs.AreaABCDEFGHB0.30841C0.225670.01147D0.166740.142180.10766E0.209690.315890.271090.05704F0.215690.114620.088800.040660.23808G0.347550.070340.049930.255500.430620.14424H0.238640.475260.385880.163290.056060.392590.52454I0.35317*−0.008920.039300.242080.438520.16998*0.026830.56396**Significant differences after Bonferroni correction (*p < 0.001)

Fig. 3. Genetic population structure of Yaku sika deer (Cervus nippon yakushimae) and thickness of the pyroclastic flow deposited on Yakushima. Black triangle (▲) represents the highest elevation point (1,936 m asl). Thickness of pyroclastic flow deposits estimated by spatial interpolation (Gaussian model, 5 km mesh) based on Geshi^36^. SAMOVA ver. 2.0^92^ was used to divide the Yaku sika deer into two genetic populations, i.e., north (open circles ○; areas B/C/G/I) and the other areas (filled circles ●; areas A/D/E/F/H: Fig. 1b). However, Area A, surrounded by the dashed line, did not belong to the northern group, even though it was heavily affected by the pyroclastic flow.

Demography of Yaku sika deer population on Yakushima Island

Table 3. Tajima’s D, Fu’s Fs, Raggedness index r, and Sum of squared deviation (SSD) for southern and northern populations based on the mitochondrial DNA control region 894 base pairs.AreaTajima’s DFu’s FsRaggedness indexSSD D p valueFsp value r p value SSD p valueNorthern: B/C/G/I0.13950.634−0.49360.4710.02420.8700.00820.690Southern: A/D/E/F/H0.25910.661−1.47670.3630.07070.0100.01920.010Whole0.66850.783−1.67020.3380.02590.2500.00410.100

The values of Tajima’s D and Fu’s Fs from 1,000 repeats with coalescent simulations were 0.6685 (p = 0.783) and −1.6702 (p = 0.338) for the whole population, 0.1395 (p = 0.634) and −0.4936 (p = 0.471) for the northern population, and 0.2591 (p = 0.661) and −1.4767 (p = 0.363) for the southern population, with none reaching statistical significance (Table 3). Mismatch analysis showed that the whole population and southern population had a unimodal distribution with a peak at three substitutions (Fig. 4 left and right), while the northern population showed a bimodal distribution with peak at one and four substitutions (Fig. 4 center). The sum of squared deviations (SSD) statistic and raggedness index for the whole population and northern population indicated that the observed frequency distribution fit the distribution expected under the demographic expansion model (SSD = 0.0041 and 0.0082, p = 0.100 and 0.690, raggedness index r = 0.0259 and 0.0242, p = 0.250 and 0.870; Table 3); however, it was rejected for the southern population (SSD = 0.0192, p = 0.010, raggedness index r = 0.0707, p = 0.010; Table 3).

Fig. 4. Mismatch distribution of Yaku sika deer (Cervus nippon yakushimae) based on the mitochondrial DNA control region 894 base pairs. The thick solid line represents the observed mismatch distribution (—●—). The thin solid line represents the demographic expansion model (—□— ). Two dashed lines indicate 95% confidence intervals (----).

Haplotype network and molecular phylogenetic tree of Japanese sika deer

The haplotype network diagram drawn using the Yaku sika deer sequences obtained in this study and the 91 sika deer sequences (Supplementary Table S4 online) registered at the NCBI showed that the populations in Japan could be clearly divided into two lineages, north and south, which is similar to the findings of previous studies (Fig. 2). The Yakushima and Kuchinoerabujima populations formed a single cluster at the edge of the southern lineage, with the Tanegashima and Mageshima populations located nearby and the Mainland Kyushu, Tsushima/Goto, Kerama, Shikoku, Honshu, and Hokkaido (Fig. 1a,b) populations located more distally (Fig. 2). All populations west of Mainland Kyushu were included in the southern lineage, while the Shikoku and Honshu populations were divided into both northern and southern lineages.

The molecular phylogenetic tree (Fig. 5) showed that after the sika deer split into the north-south lineage, Yaku sika deer (inhabiting Yakushima and Kuchinoerabujima) and Mage sika deer (inhabiting Tanegashima and Mageshima) diverged from the south lineage and further diverged into the current Yaku sika deer and Mage sika deer.

Based on divergence time estimates of the north–south lineage by Nagata et al.^12^, Kuwayama and Ozawa^19^, and Liu et al.^20^ —0.35 mya, 0.41 mya, and 0.53 mya, respectively—the divergence of this lineage was calibrated at 0.43 mya. Under this calibration, the Yakushima populations diverged from the Honshu population at 0.114 mya, and subsequently, Tanegashima and Yakushima themselves split approximately 37,500 years ago. A simultaneous radiation within the Yakushima population is estimated to have occurred sometime after 4,900–6,500 years ago.

Fig. 5. Maximum likelihood tree based on the mitochondrial DNA control region for the 110 haplotypes of the sika deer (Cervus nippon) in Japan with Wapiti (C. canadensis) as the outgroup. Numbers of nucleotide substitutions per site indicated by scale are the Tamura 3-parameter distances^97^. Numbers on branches indicate bootstrap values after 1,000 replicates. See Supplementary Table S4 for haplotype information.

Discussion

Analysis of a partial control region mtDNA sequence of Yaku sika deer revealed that the haplotype network shows a radial shape centered on Yk01^41^, with most haplotypes connected by a single mutation (Fig. 2a). The mismatch distributions of the whole Yakushima population were unimodal with a peak at three (Fig. 4 left), suggesting recent simultaneous dispersal of subspecies of sika deer on Yakushima^42,43^. The SSD and ruggedness index supported this finding. Of the four categories classified by Grant and Bowen^44^ based on different combinations of small and large haplotype diversity and nucleotide diversity values, populations that showed high haplotype diversity (h > 0.5) and low nucleotide diversity (π < 0.005) indicated that expansion occurred after a sustained period of small population size. The sika deer population in Yakushima met this criterion (Table 1). Moreover, such a demographic history was suggested by the results of the haplotype network and mismatch distribution. However, Tajima’s D (p = 0.783) and Fu’s Fs (p = 0.338) results did not suggest significant population expansion, contraction, or subdivision (Table 3). These results can be interpreted as the Yaku sika deer population sufficiently recovering from past population reductions and subsequently undergoing further generational changes that reduced the bottleneck and founder effects. In addition, unevenly distributed haplotypes on the island indicate that distinct genetic differentiation occurred among the areas (Table 2, Supplementary Tables S2 and S3 online). Individual Yaku sika deer have a much smaller home range size and are more sedentary than other subspecies of sika deer^33,35^. The presence of predators usually induces prey species to disperse widely and migrate frequently to avoid predation^45–47^; however, this was not observed for the Yaku sika deer in Yakushima, where predators have been absent in the past. Thus, the low genetic exchange between various areas in Yakushima might have induced spatial differences in haplotype structure after population reductions.

SAMOVA revealed two genetic structural populations: northern (Areas B, C, G, and I) and southern (Areas A, D, E, F, and H) (Fig. 3). The boundary of the populations roughly coincides with the thickness of the sediments from the Koya pyroclastic flow, suggesting that the extreme reduction in the Yaku sika deer population (bottleneck) was caused by the Koya pyroclastic flow from the eruption of the Kikai caldera 7,300 years ago^36^. However, the population in Area A did not belong to the northern population, despite Area A being severely affected by the pyroclastic flow. Terada et al.^17^ analyzed microsatellite DNA using a spatially explicit Bayesian model and divided Yaku sika deer into four populations: northern, southern, eastern, and western (which approximately corresponded to Area A). The population in Area A was genetically different from the populations in any other area (Table 2). Therefore, the classification of this population as an independent genetic population is reasonable. People have settled in the few flatlands in low-elevation areas around the steep island of Yakushima. In contrast, permanent settlements have not formed in Area A, which is likely because the characteristics of Area A, including low precipitation^48^, low soil moisture^49^, and particularly low temperatures among the lower elevation areas of the island^50^. Area A is the only lowland designated as a World Heritage Area based on the continuous forest extending from the mountain tops to the coast, which maintains an extremely high deer population density. These conditions may have created a unique habitat that influenced the genetic differentiation of Area A. Genetic diversity in the north was high despite the high impact of the pyroclastic flow; conversely, diversity was relatively low in the south (Table 1). This finding indicates that the ancestral haplotype Yk01 was maintained in a stable subpopulation for a long time in the southern areas, where damage from the pyroclastic flow was low. The SSD and ruggedness index also suggest that the southern population did not experience a bottleneck (Table 3). In contrast, the northern areas were severely damaged by the pyroclastic flow, although very small populations could have survived separately in the intricate mountainous terrain.

Mismatch analysis showed that the northern population had a bimodal distribution with a peak at one and four substitutions (Fig. 4 center). Bimodal mismatch distribution is a common pattern among secondary contacts in fragmented populations^51^. However, given the proximity of the two peaks in the northern population, it is more likely that multiple haplotypes remained after the bottleneck event rather than a possible secondary contact. Indeed, the SSD and raggedness index also indicate that the observed frequency distribution fits the distribution expected under the demographic expansion model. The founder population may have grown faster than the subsequent recovery population^6^. In ungulate species, early recovery of genetic diversity has been observed in founder populations^6,52^. Moreover, genetic diversity in founder populations is higher than that in populations that did not experience recent population reductions^53^. Thus, the small populations in the northern areas (founder populations) may have recovered rapidly with vegetation recovery. Kimura et al.^54^ conducted a radiocarbon analysis of pollen from buried soils at approximately 1,000 m above sea level (asl) in Area H, where the impact of the pyroclastic flow was relatively small, and the results suggested that grassland communities dominated by grasses and ferns were established at some point during the 1,100 years following the pyroclastic flow. Even in the southernmost part of Kyushu (Satsuma and Ohsumi Peninsulas, located approximately 60 km from Yakushima), which was affected by the Koya pyroclastic flow, the previously distributed evergreen broadleaf forests of Castanopsis and Lauraceae trees disappeared and did not recover for 600–900 years^55^. Thus, these findings demonstrate that the forest vegetation on Yakushima was almost completely destroyed except in certain south areas of the island and that vast grasslands or sparse forests spread for approximately 1,000 years afterward. With vegetation recovery, phytophagous animal populations should have increased based on the abundant food resources^6,53,56^. The small population that survived in the northern part of the island probably expanded rapidly as the vegetation recovered, and radiated following approximately 6,000 years ago. This finding is almost consistent with the estimated timing of the simultaneous radiation of Yaku sika deer (approximately 4,900–6,500 years ago), which was calculated by calibrating the results of the molecular phylogenetic analysis (Fig. 5) using the divergence date of 0.43 mya between the northern and southern lineages. Yaku sika deer may have acquired high genetic diversity through generational changes and genetic drift sufficient to fix region-specific mutations.

Studies on the strait formation history of these regions suggest that ancient Yakushima—that is, Yakushima prior to its separation from Tanegashima and Magejima—separated from the Paleo-Honshu Island (present-day Honshu, Shikoku, and Mainland Kyushu) 0.1–0.15 mya around MIS5^57,58^. Subsequently, some researchers propose that ancient Yakushima reconnected with Honshu during the Last Glacial Maximum (MIS2)^59,60^, while another posits that the islands remained unconnected following strait formation at 0.1–0.15 mya^57^. Molecular phylogenetic analysis conducted in this study (Fig. 5) suggests that ancient Yakushima was not connected to Honshu during the last glacial period, because of the estimated divergence of Yaku sika deer and Mage sika deer from other southern Japanese sika deer lineages at approximately 0.114 mya.

Although the Yakushima population belong to the southern Japanese lineage in the mtDNA control region, their karyotype is not the same as other populations in the southern lineage (2n = 64 or 66 chromosomes) but consistent with that of the northern lineage (2n = 68)^21^. In karyotypic polymorphisms of chromosomes in Cervidae, 2n = 68 is the basic type^21,61^ and is ancestral to 2n = 64 and 2n = 66.

The mtDNA control region of sika deer contains repetitive sequences consisting of 38–40 bp called variable number tandem repeated sequences (VNTRs)^12,62,63^. The VNTR copy numbers observed in the sika deer analyzed in this study ranged from three to ten variants, including the original sequence. (Supplementary Table S4 online). The number of VNTR variations in mammals occurs through replication slippages of the original sequence^12,62,64^, with lower repeat copy numbers representing a more ancestral form. Sika deer likely arrived in the Japanese Archipelago from the Eurasian Continent across a land bridge over the East China Sea around 0.43 mya (MIS12)^65^, which is close to the time when sika deer are thought to have split into the north-south lineage^12,19,20^. At that time, both lineages likely had 2n = 68 karyotypes and three VNTRs. After ancient Yakushima Islands separated from the Paleo-Honshu Island (present-day Honshu + Shikoku + Mainland Kyushu) at 0.1–0.15 mya^57,58^, a new southern lineage population with 2n = 64 and 66 chromosomes may have been originated in the Paleo-Honshu Island (present-day Honshu + Shikoku + Mainland Kyushu). In addition, the old type southern lineage population with 2n = 68 chromosomes was likely replaced by the new southern lineage population but remained only on the island of ancient Yakushima.

The origin of the deer inhabiting Kuchinoerabujima is a mystery. Kuchinoerabujima Island is estimated to have formed at approximately 0.58 mya based on dating the oldest tephra found on the island, the Koseda pyroclastic flow, using the glass isothermal plateau fission track method^66,67^. Kuchinoerabujima is surrounded by a sea with a depth of 500–600 m^66^, which indicates that the island has never been connected to any other island or continent since its formation. Of the two haplotypes obtained from deer on Kuchinoerabujima in this study, Kc01 was identical to KUE201–KUE206 (Supplementary Table S4 online) reported by Nagata et al.^68^ (in an 894 bp comparison). Although the read length of Kc02 was rather short (841 bp), the sequence was identical to Yk06 of Yakushima, and the haplotype was detected for the first time on Kuchinoerabujima. Sika deer can swim in the sea and have been reported to cross about 3 km from island to island^69,70^. The distance between Yakushima and Kuchinoerabujima is approximately 12 km. Although the possibility of deer swimming between the islands cannot be excluded, it is unlikely because much more frequent migration should have occurred between Yakushima and Tanegashima. The sea depth between Yakushima and Tanegashima is less than 130 m^57^, and when the sea level dropped during MIS2 (estimated from the bathymetry in Inouchi^57^, the distance between islands likely reached 2.5 km (80 m depth drop) to 5 km (60 m depth drop). However, the haplotype network (Fig. 2) and molecular phylogenetic tree (Fig. 5) clearly distinguish populations of Yaku and Mage sika deer, and common haplotypes have not been found between the two islands. (see also Nagata^25^ and Yamada et al.^14^). Therefore, it is unlikely that sika deer (at least female) migrated between Yakushima and Tanegashima. Calibration of the molecular phylogenetic analysis (Fig. 5) indicated that the Kc01 haplotype, found exclusively on Kuchinoerabujima, was formed approximately 2,100 years ago, during the Yayoi period; therefore, it is early to infer whether they were introduced by humans after Antiquity. Further detailed genetic analyses, including nuclear DNA analyses, are required to address the origin of the deer at Kuchinoerabujima.

In this study, we demonstrated the recovery process of genetic diversity and the reconstruction of local genetic structures after a bottleneck within an island population of deer. Different degrees of past catastrophic impact areas across the island reflected the current local genetic diversity and the mosaic composition of the population structure. These findings suggest these genetically distinct subpopulations—particularly the Area A population—can be treated as independent conservation management units, which is essential for the development of effective long-term monitoring and conservation strategies. The process identified in this study might also have occurred in other populations, thus leading to an early recovery of genetic diversity.

Materials and methods

Study sites

Yakushima (30° N, 130° E) is a round island with an area of 504 km^2^ and a circumference of 130 km. It is located 60 km south of Kyushu Island, southern Japan, and lies 18 km west of Tanegashima Island (Fig. 1a,b). Most of the island is steep and mountainous, with Mt. Miyanoura (1,936 m asl) standing at the center of the island and representing the highest mountain in the Kyushu area. The average annual precipitation in Yakushima shows great regional variation, ranging from 2,400 mm in the western coastal region (lowest) to 7,400 mm in the eastern inland region (highest)^48^. The western part of the island presents a continuous vertical distribution of natural vegetation from subtropical evergreen broad-leaved forests in coastal areas to cool-temperate grasslands at the summit, and its unique ecosystem was designated as a World Natural Heritage Site (10,747 ha) in 1993. More than 90% of Yakushima is covered with forest^71^, and most agricultural, residential, and commercial areas are located in low-lying lands at an altitude of less than 150 m asl^72^. Kuchinoerabujima (36 km^2^) is a volcanic island located 12 km west of Yakushima Island. When Mt. Shindake (626 m asl) erupted in 2015, approximately 150 islanders were forced to evacuate from the island, and their return was only possible after seven months. Today, approximately 40% of Yakushima and the entire area of Kuchinoerabujima together form Yakushima National Park.

Between 1.0 and 1.7 mya, Yakushima, along with Tanegashima, was likely connected to Paleo-Honshu Island (present-day Honshu + Shikoku and Mainland Kyushu; Fig. 1a) and the Eurasian Continent^59,73^. According to marine sediment surveys concerning the formation of the strait, an island detached from Honshu—which comprises present-day Shikoku and Kyushu—between 0.15 and 0.1 mya^57,58^, and was subsequently divided into Tanegashima and Yakushima (Fig. 1b). Many remains of the early Jomon period have been found in Yakushima, where human settlements are estimated to have formed at least 7,000 years ago^74^. At approximately 0.58 mya, the volcanic island Kuchinoerabujima was created^66,67^, and was likely never connected to the surrounding islands or the continent (Fig. 1b). The Kikai Caldera, located approximately 30 km northwest of Yakushima (Fig. 1a,b), erupted approximately 7,300 years ago, and a large amount of pyroclastic sediment covered Yakushima and Kuchinoerabujima^36^.

Sample collection, storage, and DNA extraction

We mainly used deer feces as noninvasive samples for genetic analysis. Fecal samples are the most informative and highly available sources of Yaku sika deer genomic DNA^75^. This approach enabled the collection of genetic material even from steep mountainous regions and wildlife protection areas without relying on conventional invasive methods (e.g., tissue samples from harvested deer). Fresh fecal samples are desirable for obtaining accurate genetic information; however, the large regional differences in deer densities on Yakushima^76,77^ increase the difficulty of collecting fresh feces in low-density areas. However, analyses of brown bear (Ursus arctos) feces collected more than one month after defecation^78^ and Asian elephants (Elephas maximus) collected eight days after defecation^79^ revealed that the mtDNA levels were sufficient for analysis. Our previous study on Yaku sika deer showed that feces collected within three days after defecation can be used to analyze nuclear microsatellite DNA regardless of temperature or rainfall^80^. Because mtDNA has overwhelmingly more copies per cell than nuclear DNA, it generally retains its analytical efficiency longer than nuclear microsatellite DNA^81^. Thus, feces collected several days after defecation was acceptable for our analysis.

Genomic DNA was collected from Yaku sika deer from February 2004 to April 2017. We collected 373 fecal samples—in addition to 32 hair samples and 14 tissue fragments—from Yakushima. We also collected five fecal samples from Kuchinoerabujima in July 2004. In the western lowlands of Yakushima, we performed deer tracking and collected fresh feces immediately after defecation because deer in this area are habituated to humans and are not afraid of researchers; therefore, direct observations are relatively easy^31,32^. In other areas, we walked through the forest to search for fecal remains on the ground and collected only feces that appeared to be relatively fresh based on their appearance, as suggested by Agetsuma-Yanagihara et al.^80^. In addition, we used hairs collected from deer captured for transmitter installation^31,32,35^, fallen hair bundles in the forest, and tissue fragments collected from carcasses found in the forest and from individuals culled as pests by the local government. Deer capture and handling protocols for the transmitter installation had been reviewed and permitted annually by the Ministry of the Environment, Government of Japan and the Kagoshima Prefectural Government throughout the study period.

Feces collected from February 2004 to May 2013 were air-dried at 15 °C to 25 °C and stored at −20 °C to −30 °C. The surface layer of the frozen feces was scraped with a sterile scalpel on ice, and DNA was extracted and purified using the QIAamp DNA Stool Mini Kit (QIAGEN, Hilden, Germany). From May 2015 to April 2017, the surface of feces, which is thought to contain adhered intestinal cells^82^, was rubbed several times with a sterile cotton swab, which was then rinsed in a 2 mL tube containing conditioned lysis buffer^83^. The solution was stored at room temperature in a cool, dark place^80^. DNA was extracted and purified from the fecal solutions using the QIAamp DNA Stool Mini Kit (QIAGEN) or the QIAamp Fast DNA Stool Mini Kit (QIAGEN). Hair samples were washed thoroughly with ethanol (> 70%) and stored at −20 °C to −30 °C, and DNA was extracted and purified using ISOHAIR (Nippon Gene Co., Ltd., Tokyo, Japan). Immersed tissue fragments were stored in ethanol (> 70%) at −20 °C to −30 °C, and DNA was extracted and purified using a DNeasy Blood & Tissue Kit (QIAGEN). The purified DNA was stored at −20 °C to −30 °C if not analyzed immediately.

PCR direct sequencing of the mtDNA control region

The control region of mtDNA was amplified using the primer pair LD5 (5′-AAGCCATAGCCCCACTATCAA-3′) and H597 (5′-AGGCATTTTCAGTGCCTTGCTTTG-3′)^84^ and the Tks Gflex DNA Polymerase PCR enzyme (Takara Bio Inc., Shiga, Japan), which is an alpha-type enzyme with both proofreading and polymerase activity^85,86^. In addition, previous study revealed that Tks Gflex has high fidelity^87^. PCR amplifications were carried out in a total volume of 20 µL, which contained 2 µL template DNA, 0.4 µL Tks Gflex DNA Polymerase, and 0.05 µL of each 100 µM primer (LD5 and H597). PCR was performed on a thermal cycler (TaKaRa PCR Thermal Cycler Dice Touch; Takara Bio Inc.). The following PCR protocol was used: after initial incubation at 94 °C for 1 min, 45 amplification cycles were performed with denaturation at 98 °C for 10 s, annealing at 55 °C for 15 s, and extension at 68 °C for 40 s. Successful PCR amplification was confirmed by electrophoresis of PCR products on a 1.5% agarose gel containing GelGreen (Biotium, Fremont, CA, USA), followed by visualization under UV illumination.

PCR products were purified with the High Pure PCR Product Purification Kit (Roche Diagnostics K.K., Tokyo, Japan) and sequenced using the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Waltham, MA, USA). Cycle sequencing reactions were carried out in a total volume of 10 µL, which contained 2 µL PCR product, 1 µL Ready Reaction Mix, 2 µL Sequencing Buffer, and 1 µL primer (1 µM; LD5 or H597). After initial incubation at 96 °C for 1 min, 50 amplification cycles were performed with denaturation at 96 °C for 10 s, annealing at 50 °C for 5 s, and extension at 60 °C for 2 min. DNA fragments were sequenced using an ABI3130 Genetic Analyzer (Applied Biosystems).

Data analysis

We visually confirmed the sequencing chromatograms using MEGAX v10.1.6^88^ and determined the nucleotide sequence of each sample. Haplotypes for each sample were determined using DnaSP ver. 6.12.03^89^. Fecal samples collected within a 200 m radius that exhibited completely identical sequences and could not be determined to originate from different individuals based on external appearance, such as size and shape, were excluded from subsequent genetic analyses. For the haplotypes obtained in this study (18 from Yakushima and 1 from Kuchinoerabujima excluding Kc02) and those found in the sequence data of sika deer in Japan registered in the NCBI (Supplementary Table S4 online), we reconstructed median-joining networks^90^ using Network 10.2.0.0 (https://www.fluxus-engineering.com/sharenet.htm). The positional definition and counting method of VNTRs differ among studies by Randi et al.^63^, Nagata et al.^12^, and Cook et al.^62^; in the present study, we followed the sites and method of Randi et al.^63^. A total of 42 distinct VNTR units were identified from all sequence data used in the network analysis, with repeat copy numbers (including the original) ranging from three to ten. The original VNTR sequence and the two copies are commonly found across all subspecies of sika deer, as well as in the Manchurian and North American subspecies of elk (C. canadensis), rusa deer (C. timorensis), and sambar deer (C. unicolor)^63^. In contrast, the presence of three or more copies is a unique feature observed only in sika deer populations in the Japanese Archipelago^12,63^. Thus, the original and the two-copy sequences represent ancestral forms, and subsequent copies are considered to be derived variants. Accordingly, we used only the original VNTR sequence and its two copies in the haplotype network and molecular phylogenetic analyses in this study. The sites where the bases could not be deciphered were excluded from the analysis, and the alignment gap was treated as a fifth base.

To examine the genetic subpopulation structure of Yaku sika deer in Yakushima, we divided Yakushima into nine areas from A to I (Fig. 1b) based on the following habitat conditions: impact of the pyroclastic flow caused by the giant eruption of the Kikai Caldera approximately 7,300 years ago^36^; forest logging history over the past 100 years^91^; distribution of annual precipitation^48^; altitude; and status as a World Heritage Area (strictest natural protection area). We also searched for genetic boundaries in the aforementioned nine areas based on SAMOVA ver. 2.0^92^. To examine the relationship between the results of SAMOVA and impact of the Koya pyroclastic flow, the distribution of the sediment thickness over Yakushima was estimated by spatially interpolating (Gaussian model, 5 km mesh) the thickness measured by Geshi^36^ using ArcGIS ver. 9.3 (ESRI Inc.).

We employed Arlequin v3.5.2^93^ to calculate the haplotype frequency, haplotype diversity (h), nucleotide diversity (π), and pairwise FST values for each area and performed AMOVA. We computed Tajima’s D^94^ and Fu’s Fs^95^ values to summarize the haplotype frequency, and their randomness and neutrality were tested using Arlequin v3.5.2. Mismatch analysis^42,43^ was performed using Arlequin v3.5.2^93^ to examine how each sample population formed the current population. The SSD between the observed and expected distributions and ruggedness index (r)^96^ were used to assess the goodness of fit of the sudden expansion model in the mismatch analysis. Moreover, we used MEGA X v10.1.6^88^ to construct phylogenetic trees of sika deer in Japan (Supplementary Table S4 online) with Wapiti (C. canadensis) as the outgroup. The best nucleotide substitution model for phylogenetic tree construction was estimated using the maximum likelihood method, which was complemented by MEGA X v10.1.6^88^.

Supplementary Information

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2