Sample and data collection

We analyzed the entire Cytb coding region (1140 bp) and CR hypervariable region I (552 bp after exclusion of indels) sequences of 106 A. speciosus individuals. Sequences were obtained from the DNA Data Bank of Japan (DDBJ), European Nucleotide Archive (ENA), and GenBank nucleotide databases (Suzuki et al. 2015); other sequences were newly determined for 43 individuals collected in Hokkaido in 2017 (Table 1). Field collection was performed with permission from the Hokkaido prefectural government for academic research only, following the Guidelines for the Procedure of Obtaining Mammal Specimens as Approved by the Mammal Society of Japan.

Table 1.

List of specimen of Apodemus speciosus used in this study

Sequencing analysis

Genomic DNA was extracted from ethanol-preserved tissue samples using a QIAamp DNA Mini Kit (QIAGEN) following the manufacturer's protocol. Cytb and CR sequences were determined as reported previously (Suzuki et al. 2015). The sequences were aligned using the Muscle function of the MEGA7 software (Kumar et al. 2016). All sequences obtained in this study were deposited in the International DDBJ/ENA/GenBank DNA database (accession nos. LC655789–LC655874).

Phylogenetic and population genetic analysis

A maximum-likelihood (ML) tree was constructed with concatenated mtDNA sequences of Cytb and CR using MEGA7 (Kumar et al. 2016) with the GTR + G + I model. The best fit model was determined using Akaike Information Criterion (AIC) implemented in MEGA7. Node reliability was assessed using 500 bootstrap replicates. Median-joining (MJ) networks were constructed using the concatenated mtDNA sequences with PopART v1.7 (Leigh and Bryant 2015).

Geographical trends in genetic diversity were examined using landscape shape interpolation analysis, implemented in Alleles In Space (AIS; Miller 2005). We converted the latitude and longitude of each sampling locality to XY coordinates on an arbitrary scale, using the plane rectangular coordinate system of Japan (Zone 12, 44°00′N, 142°25′E), following the instruction of AIS. Because sample sizes differed among localities, we used one sequence per haplotype at each location to avoid sampling bias. Using the geographic coordinates and DNA sequences of Cytb, we generated a Delaunay triangulation-based connectivity network of sampling localities and their relative genetic distances, which decreased with increasing geographic distance.

We used ARLEQUIN v3.5.1 to calculate nucleotide (π) and haplotype diversity (Hd) using each locus data (Excoffier and Lischer 2010). Mismatch distribution analyses and neutrality tests, Tajima's D, and Fu's Fs, were performed to evaluate rapid population expansion using ARLEQUIN. The significance of neutrality was assessed using 1000 replicates of coalescent simulations. Historical population expansion was assessed through mismatch distribution analysis (Rogers and Harpending 1992). The expected distribution was simulated using the sudden expansion model (Schneider and Excoffier 1999; Excoffier 2004), which was validated using a parametric bootstrap approach with 1000 replicates. For each replicate, we compared the sum of squared deviations (SSD) between the observed and expected distributions with the SSD between the simulated and expected distributions using ARLEQUIN. We applied the raggedness index (r) (Harpending 1994) as a test statistic for sudden expansion model estimation. The expansion parameter tau (τ) was estimated in each set of haplotypes (we call it Group) to search for signs of sudden demographic expansion using ARLEQUIN (Excoffier and Lischer 2010). The start times (t) of expansion events were calculated using the equation t = τ/2µk, where µ is the evolutionary rate in years and k is the length of the analyzed fragment (1140 bp) (Suzuki et al. 2015; Hanazaki et al. 2017). There is growing evidence to support that the rate of evolution of mtDNA in animals is not constant and is time-dependent, especially in the early stages of divergence, partly due to the involvement of mild deleterious mutations (e.g., Ho et al. 2005, 2011; Inoue and Suzuki 2022). Hence, here we used the time-dependent evolutionary rate of Cytb in small mammals (Honda et al. 2019; Suzuki 2021).

Phylogeny and spatial distribution of mtDNA variation

We determined the Cytb (1140 bp) and CR (552 bp) sequences of mtDNA in 43 samples of A. speciosus newly collected from Hokkaido. A 2-bp insertion observed in two CR sequences (HS5444, HS378) was removed from subsequent phylogenetic analyses. We combined these sequences with those from databases (n = 63) to obtain 106 sequences in total: 83 sequences from 23 localities on Hokkaido and 23 sequences from 20 localities on Honshu, Shikoku, and Kyushu (Table 1, Fig. 1A).

We constructed an ML tree (Fig. 2) and MJ network (Fig. 3) using the concatenated sequences of Cytb and CR. In the MJ network, we found star-shaped clusters representing Hokkaido and Honshu/Shikoku/Kyushu, as previously reported (Suzuki et al. 2015). The Hokkaido cluster encompassed a cluster with small numbers of nucleotide substitutions (Group A; Fig. 3). The remaining Hokkaido haplotypes and the Honshu cluster were designated Groups B and C, respectively (Fig. 3). The ML tree showed that neither of monophyly for Group A and B from Hokkaido was supported (Fig. 2). The geographic distribution of Groups A and B is shown in Fig. 1A. Group A was observed throughout Hokkaido, but was more abundant in central Hokkaido, whereas Group B appeared frequently in southern and eastern Hokkaido. The ML tree showed four subgroups (B1–B5) of haplotypes connected by short branches, supported by moderate to high bootstrap values (73–99%). Haplotypes from B1 and B2 tended to be confined to eastern Hokkaido, whereas those from B3, B4, and B5 were distributed in both southern and eastern Hokkaido. Overall, members of Group B, which was derived from an older period than those of Group A, tended to remain within eastern Hokkaido.

Fig. 1.

(A) Sampling localities for the large Japanese wood mouse Apodemus speciosus in Hokkaido, the northern island of the Japanese archipelago. Geographic distribution of mitochondrial Groups A and B. Dotted circle indicates the approximate area of central Hokkaido. The geographic ranges of subgroups B1–5 of Group B are shown. Dark gray shading indicates the Hidaka mountain range. (B) Genetic landscape interpolation analysis results for mitochondrial Groups A and B. Light (yellow in online version) and dark (blue in online version) shading indicate high and low genetic diversity, respectively.

We performed landscape shape interpolation analysis of the Cytb sequences from Hokkaido (Fig. 1B). In Group A, overall high genetic diversity in Hokkaido was observed, while only the southwestern part was colored dark gray (blue in online version), showing low genetic diversity. Group B showed relatively high genetic diversity in the eastern part.

Analyses for rapid population expansion

In each gene sequence dataset of Cytb and CR, the values of Hd, π, and τ were lower and higher in Groups A and C, respectively (Table 2). The value of τ was lower in Group B than Group C in the Cytb sequences, but nearly the same in Groups B and C in the CR sequences. Tajima's D and Fu's Fs neutrality tests were significantly negative in almost all datasets examined, but Tajima's D was not significant in the CR dataset of Group B. Mismatch distribution analysis showed that the sudden expansion model was not rejected in all datasets except Cytb sequences in Group B in terms of both SSD and RI (Table 2, Fig. 4).

It has been reported that the different evolutionary rate should be applied to time estimations at different ages because of the time dependency of the evolutionary rates (Suzuki 2021). Times of the onset of expansion events were therefore estimated using the Cytb data, considering the time-dependent evolutionary rate curves as shown in Cytb of small mammals, which are similar among different lineages, such as wood mice, voles, and moles (Suzuki 2021). Accounting the degree of the τ values for Cytb in Groups A (τ = 2.88), B (τ = 4.99), and C (τ = 7.84) (Table 2), we selected different evolutionary rates of 0.11, 0.047, and 0.028 substitutions/site/million years, respectively, according to the evolutionary rates estimated for each time scale of the τ value (Honda et al. 2019; Suzuki 2021). The start times of the expansion events in Groups A, B, and C were estimated to be 11 500, 46 600, and 122 000 years ago, corresponding to the early stages of MIS 1, MIS 3, and MIS 5, respectively. Furthermore, when the data set (n = 83) of Groups A plus B was examined, signals of rapid population expansion were obtained in the test of Tajima's D and Fu's Fs, and the value of τ obtained (τ = 4.41) was similar to that of Group B (Table 2), suggesting that the expansion started at the early MIS 3.

Fig. 2.

Maximum-likelihood trees constructed from 106 concatenated sequences (cytochrome b and the control region, 1692 bp) of Apodemus speciosus. Three clusters represented by haplotypes from Hokkaido (Groups A and B) and Honshu/Shikoku/Kyushu (Group C), and four subgroups of Group B, are indicated. In Groups A and B, numbers in parentheses next to species names indicate the sampling locality (Table 1, Fig. 1A). Node numbers indicate bootstrap values (1000 replicates). Bar indicates substitution numbers per site.

Possible factors shaping rapid expansion events in small terrestrial mammals

In this study, we constructed a phylogenetic tree using the Cytb and CR sequences (1692 bp) of A. speciosus from the Japanese archipelago, revealing three groups (A–C) as observed in our previous study (Fig. 2; Suzuki et al. 2015). The constructed network of the concatenated (Fig. 3) and Cytb (1140 bp; data not shown) sequences exhibited small (Group A: τ = 2.88 in Cytb), medium (Group B: τ = 4.99 in Cytb), and large (Group C: τ = 7.84 in Cytb) star-shaped characteristics in the network. The rapid expansion of Groups A and C was supported by non-significant SSD and RI statistics and significant P values for Tajima's D and Fu's Fs tests in both Cytb and CR (Table 2, Fig. 4). The rapid population expansion of Group B was generally supported, but not by SSD or RI analyses of the Cytb data or Tajima's D for the CR data (Table 2). This significant divergence between observed and expected mismatch distributions may be explained by a secondary expansion event, as indicated by subgroups B1–5 (Fig. 3).

The formation of these small and large groups was likely caused by rapid expansion events initiated 10 000 and 130 000 years ago, as post-LGM and post-PGM events, respectively (Suzuki et al. 2015). The substantial effects of these glacial maxima on population expansion have been suggested in mtDNA analyses of several Japanese mammals including the small Japanese wood mouse (Suzuki et al. 2015; Hanazaki et al. 2017), western Japanese mole Mogera wogura (Nakamoto et al. 2021), Japanese giant flying squirrel Petaurista leucogenys (Oshida et al. 2009), and red-backed vole species (genus Myodes, Honda et al. 2019). Therefore, we can conclude that the effects of both glacial maxima had a great impact on the population dynamics of A. speciosus and other small mammals in the Japanese archipelago.

Fig. 3.

Median-joining networks constructed from 106 concatenated sequences (cytochrome b and the control region, 1692 bp) of Apodemus speciosus. Three clusters are represented by haplotypes from Hokkaido (Groups A and B) and Honshu/Shikoku/Kyushu (Group C). Two subgroups (A1, A2) of Group A and five subgroups (B1–5) of Group B are indicated. Mutation numbers between haplotypes are indicated by hatch marks. Circle size is proportional to sample size.

Table 2.

Genetic diversity indices and detection of rapid expansion events with the mismatch distribution analysis of two mitochondrial DNA markers (Cytb and CR) for Apodemus speciosus

Fig. 4.

Mismatch distribution of the mitochondrial gene sequences: cytochrome b (Cytb; A–C) and the control region (CR; D–F). Each Group of A (A, D), B (B, E), and C (C, F) was analyzed. Bars indicate observed frequency; line indicates the expected frequency under the sudden expansion model. SSD, sum of squared deviations; RI, Harpending's raggedness index.

The τ values of the groups (Table 2) suggest that the medium-sized Group B (Cytb, τ = 4.99) is the consequence of an event that occurred between the post-LGM (Group A; Cytb, τ = 2.88) and post-PGM (Group C; Cytb, τ = 7.84) events. The climate change from MIS 4 to MIS 3 is known to have had a tremendous impact on the global climate (Doughty et al. 2021; Geibert et al. 2021). Rapid population growth of the Holarctic northern red-backed vole (Myodes rutilus) in North America is reported to have been associated with an expansion event ca. 50 000 years ago, when large ice sheets shrank in size throughout the continent (Kohli et al. 2015). The grassland-dwelling striped field mouse A. agrarius occurring in a broad area of northern Eurasia from Europe to Asia showed a rapid expansion signal during MIS 4-3 in mtDNA sequence data (Kozyra et al. 2021). Interestingly, the forest-dwelling Asian house rat (Rattus rattus complex mitochondrial Lineage II, defined by Pagès et al. 2010 and Aplin et al. 2011), which occurs in East and Southeast Asia including Myanmar, also showed rapid expansion signals in mtDNA analysis (Maung Maung Theint et al. 2020), suggesting that the effects of climate change on small rodents associated with the early MIS 3 were widespread and numerous at the global level. In the Japanese archipelago, MIS 4 was characterized by more snowfall and greater glacial development due to the presence of high mountain ranges (e.g., Hidaka Mountains, Hokkaido; Fig. 1) compared to the climate change impacts of the LGM (Sawagaki et al. 2004), and therefore it might be sufficient to trigger a population bottleneck as a prerequisite event for an expansion event to occur (Hanazaki et al. 2017). In an mtDNA study of A. argenteus, Hanazaki et al. (2017) concluded that an expansion event was triggered by warming during the early MIS 3. In our previous mtDNA study of grassland Myodes voles, which occur in a broad area of Hokkaido (M. rufocanus) and in the high elevation grasslands of western Honshu (M. smithii), we observed signals of multiple expansions, presumably initiated during early MIS 3 (Honda et al. 2019). These considerations led to the hypothesis that changes in the environment during the late MIS 4 and early MIS 3 caused the rapid expansion of Group B in the Hokkaido population of A. speciosus.

Possible scenarios for the formation of the spatial genetic structure of A. speciosus in Hokkaido

The distribution patterns of the genetic structure of Groups A and B of A. speciosus populations in Hokkaido are shown in Figs. 1 and 5A. The newly generated Group A was found to be dominant in central Hokkaido, but coexisted with the older Group B in peripheral areas. Generally, glacial refugia are assumed to harbor more genetically distant alleles than areas colonized after a glacial period. Thus, in addition to Group A, several distinct lineages belonging to Group B were detected in the eastern part of Hokkaido (Fig. 1A), suggestive of its role as a glacial refugium. In contrast, limited numbers of distinct lineages were detected in central and southern parts of Hokkaido, suggesting that the populations of western part of Hokkaido were newly formed by westward migration from eastern Hokkaido during the post-LGM warm period.

Fig. 5.

Geographic distribution of distinct mitochondrial lineages in Japanese wood mice. Schematic representation of Groups A and B in A. speciosus from Hokkaido (A, this study), and two distinct clusters of A. speciosus from the Japanese archipelago (B, Suzuki et al. 2004; Tomozawa and Suzuki 2008; Tomozawa et al. 2014; Suzuki et al. 2015). This geographic structure is well explained by the concept that sites with low genetic diversity serve as sources of postglacial colonization, while those with high genetic diversity become sinks for migrating haplotypes. In addition, simultaneous mtDNA expansion occurs as one haplotype tends to drive out the others, whereas islands and apical regions tend to evade this effect.

Alternatively, the high genetic diversity indicated by the presence of several distinct lineages of Group B in eastern Hokkaido can also be explained by the migration of haplotypes that survived the LGM in south and central parts of Hokkaido. This can be considered in the context of source-sink dynamics (e.g., Fazio et al. 2004). It assumes an eastward shift of glacial retreat in Hokkaido, which can be supported by a pollen fossil study that showed increasing dominance of Quercus (oak trees) in southwest Hokkaido ca. 1000 years earlier than in northern Hokkaido (Igarashi et al. 2011). This indicates that it is important to note that areas of high genetic diversity are not necessarily sources of diversity, but rather may be areas that have characteristics as sinks. In A. speciosus, the mtDNA clade (Clade II) from offshore islands of Sado Island, Hokkaido Island, Izu Islands, and Satsunan Islands (Fig. 5B) are thought to have originated by intermittent migrations from the central main islands to the offshore islands during the construction of the land bridge-like structure at the periods of glacial maxima (Suzuki et al. 2004; Tomozawa and Suzuki 2008), resulting in a higher degree of genetic diversity in the offshore islands as a whole. This is explained in the context of source-sink as well.