Impact of late Quaternary climate change on the demographic history of Japanese field voles and hares revealed by mitochondrial cytochrome b sequences

Hitoshi Suzuki; Mitsuo Nunome; Takuro Yanase; Takeshi Eto; Masashi Harada; Gohta Kinoshita

doi:10.1266/ggs.24-00145

ABSTRACT

The mitochondrial cytochrome b gene (Cytb) of the Japanese field vole (Microtus montebelli), an herbivorous rodent, was subjected to an analysis of sequence variation with the objective of elucidating the population histories of this species. Construction of a phylogenetic tree revealed the existence of several region-specific lineages in Honshu and Kyushu, which were evenly separated from each other. In consideration of the documented time-dependent evolutionary rates of rodents, the estimated divergence times indicate that the region-specific lineages of M. montebelli emerged 160,000–300,000 years ago. In a haplotype network, the region-specific lineages from northern and central Honshu tended to show star-shaped clusters, with additional internal star-shaped clusters, indicative of two periods of population expansion. The onsets of these expansions were estimated to have occurred 15,000 and 10,000 years ago, respectively, suggestive of association with the two periods of rapid warming following the last glacial maximum (LGM). In contrast, such predicted post-LGM expansion events were less pronounced in the southern lineages, implying latitudinal dependence of the effect of the LGM on population dynamics. Sado Island haplotypes exhibited a network with a star-shaped pattern and a 10,000-year-old expansion signal, surrounded by a Honshu haplotype cluster with a 15,000-year-old expansion signal, suggesting that post-LGM expansion events contributed to the formation of the Sado population. A reanalysis of Cytb sequences of the Japanese hare (Lepus brachyurus), which has a similar geographic range to the voles, yielded results that were consistent with those of the vole analysis, confirming that the characteristics of the post-LGM expansion event were dependent on latitude, involved two successive expansion events, and enabled migration across deep straits. It seems reasonable to infer that the environmental changes that occurred during the warm periods following the LGM were a contributing factor in the expansion of the distribution range of newly emerged haplotype groups.

INTRODUCTION

Global climate change during the Quaternary period has significantly impacted the population dynamics of terrestrial organisms in the Northern Hemisphere, especially at high latitudes (Hewitt, 1996, 2000). This period is characterized by glacial and interglacial cycles of approximately 100,000 years (Lisiecki and Raymo, 2005). However, the full impacts of these changes on different species and populations from regions at varying latitudes are not yet fully understood. In addition, global climate change has led to fluctuations in sea levels, which have resulted in the formation and disappearance of land bridges (Lambeck et al., 2002). This has led to genetic exchange over time between regions that were once isolated by deep sea straits. Conversely, naturally occurring populations inhabit offshore islands that are separated by deep water that have never had land bridges, even during the glacial maxima. There is interest in the timing of cross-strait migration and the means by which it occurred, e.g., by land bridge formation or on floating objects (Houle, 1998; Sato, 2017; Yang et al., 2018).

Small to medium-sized mammals in the Japanese archipelago are of interest for understanding the effects of climate change on the establishment and maintenance of genetic diversity. This archipelago, which is long and narrow extending from north to south and has varying elevations, provides habitats for numerous terrestrial mammal species (e.g., Sato, 2017). Hokkaido, situated at high latitudes, harbors many subarctic species (Ohdachi et al., 1997; Kinoshita et al., 2024). In contrast, the three main mid-latitude islands—Honshu, Shikoku and Kyushu—are home to a considerable number of temperate species that persisted throughout the Quaternary period (Serizawa et al., 2000; Nunome et al., 2007; Kinoshita et al., 2019). The Japanese archipelago also includes islands separated by deep sea channels. Some straits, such as that between the Oki Islands and Honshu, are thought to have formed a land bridge during the Quaternary glacial maximum due to sea level fluctuations (Japan Association for Quaternary Research, 1987). In contrast, others, such as those between Hokkaido and Honshu and between Sado Island and Honshu, are too deep (exceeding 150 m) to have formed a land bridge even during the glacial maximum (Japan Association for Quaternary Research, 1987). Nevertheless, Sado Island is home to a number of terrestrial mammals, including the Japanese field vole (Microtus montebelli), Japanese hare (Lepus brachyurus), Sado mole (Mogera tokudae) and large Japanese wood mouse (Apodemus speciosus). Mogera tokudae and A. speciosus were shown to have distinct lineages on the island based on genetic analysis of mitochondrial DNA (mtDNA) (Tsuchiya et al., 2000; Tomozawa and Suzuki, 2008; Tomozawa et al., 2010, 2014). However, the full genetic and historical context of naturally occurring mammal species on Sado Island, including M. montebelli and L. brachyurus, remains to be elucidated.

The typical pattern of mtDNA diversity observed in regional populations is that of admixture, which can be attributed to a series of rapid population expansion events that occurred in disparate regions over a span of time (Hewitt, 2004). However, by focusing on the spatiotemporal dynamics of individual expansions or lineages, it is possible to gain a more comprehensive understanding of the population dynamic history of the species as a whole, based on a clear understanding of the evolutionary rate of mtDNA. Nevertheless, the hypothesis that the rate of evolution is contingent upon the time of divergence (Ho et al., 2011) is currently the subject of further investigation. We have investigated the evolutionary rate of the mitochondrial cytochrome b gene (Cytb) by analyzing the population dynamics of Japanese land mammals (Suzuki et al., 2015; Hanazaki et al., 2017; Honda et al., 2019; Nakamoto et al., 2021; Suzuki, 2021). In the temperate Japanese archipelago (e.g., Honshu and Kyushu), there are signals not only of post-last glacial maximum (LGM) divergence, but also of divergence that occurred at earlier times. For example, a study of mitochondrial haplotypes of Apodemus argenteus (the lesser Japanese wood mouse) and A. speciosus revealed star-shaped networks of three different sizes (large, medium and small) (Suzuki, 2021). Based on the pattern of variation in fossil pollen of oak species used as rodent resources over the past 144,000 years (Igarashi and Oba, 2006), the small, medium and large haplotype groups of A. speciosus are considered to be the result of rapid expansion events since the end of the low-temperature period of the LGM, marine isotope stage 4 (MIS 4) and the penultimate glacial maximum (PGM), respectively (Suzuki et al., 2015; Hanazaki et al., 2017). The evolutionary rates of Cytb were determined for the two Japanese Apodemus species using these periods as calibration points (Suzuki et al., 2015; Hanazaki et al., 2017). A comparable pattern of evolutionary rate in Cytb has been indicated by genetic studies of other small Japanese mammals, including voles (genus Myodes) and moles (genus Mogera) (Honda et al., 2019; Nakamoto et al., 2021; Suzuki, 2021). It can thus be inferred that this time-dependent evolutionary rate pattern is likely to be shared across various species, including hares.

Two herbivorous mammals, M. montebelli and L. brachyurus, live in the mid-latitude regions of Japan. The Microtus voles are found on Honshu, Kyushu and Sado Island (Ohdachi et al., 2015). Voles of the genus Myodes (M. andersoni in the north and M. smithii in the south) inhabit the highlands, while Microtus voles live in lowland grasslands. Although M. montebelli has attracted attention from various perspectives, such as forestry damage, ecology and zoonosis (Asakawa, 2005; Murano et al., 2023; Tsukada et al., 2023), there have been no studies on the genetic structure and population dynamics of Microtus voles. The Japanese hare L. brachyurus is found on Honshu, Shikoku and Kyushu, as well as a number of offshore islands, including Sado and Oki islands (Ohdachi et al., 2015). Previous studies identified distinct geographic lineages that show population dynamics related to Quaternary climate change (Nunome et al., 2010). However, a comprehensive analysis of spatiotemporal dynamics using the time-dependent evolutionary rates has not been performed.

The objective of this study was to examine the influence of Quaternary climate change on regional phylogenetic differentiation and population dynamics in M. montebelli. We examined intraspecific variation in the mtDNA of M. montebelli collected across a broad region in Honshu and Kyushu. To gain insight into the temporal aspects of evolutionary change, we employed a time-dependent approach to examine the evolutionary rates of rodents. Furthermore, the genetic variation in the mtDNA of voles from Sado Island was examined with the objective of gaining insights into the natural history of trans-strait vole colony formation. Furthermore, the previously determined Cytb sequence for L. brachyurus (Nunome et al., 2010), which exhibits a comparable distribution range with M. montebelli, including its presence on Sado Island, was also employed in the present study. By employing a comparative methodology, we aimed to achieve a comprehensive understanding of the impact of Quaternary climate change on the phylogeography and population dynamics of small mammals in the Japanese archipelago.

RESULTS

Cytb sequence variation in M. montebelli

We newly determined Cytb sequences (n = 50) of M. montebelli from Honshu, Kyushu and Sado Island (Fig. 1). Together with conspecific sequences from the databases (n = 3), we constructed a maximum likelihood (ML) phylogenetic tree with the TN93 + G + I model and outgroup taxon (Fig. 2A). The tree displayed a brush-like structure, with a long solitary branch extending from the outgroup lineage and multiple branches near the end of the branch. We observed several different intraspecific lineages, which we tentatively grouped into six groups called lineages I–VI. The relationships among these lineages were unresolved, and the nodes were connected by short branches, indicating their divergence within a relatively short period. Lineages II and III were characterized by a single haplotype in Niigata and Nagano prefectures, respectively, while the remaining four lineages were composed of multiple haplotypes and appeared as variants specific to certain regions (Fig. 1). In particular, the haplotypes of lineage I were prevalent in northern and central Honshu, including Sado Island. Lineage IV was identified in the Hokuriku and Kansai regions. Lineage V included haplotypes from the Kansai and Chugoku regions, while lineage VI included two distinct haplotypes found in Kyushu designated as VIa and VIb.

Fig. 1. A map of the Japanese archipelago showing the distribution range of Microtus montebelli occurring on Honshu, Kyushu and Sado Island. The geographic distribution of the six major mitochondrial cytochrome b lineages (lineages I to VI) is indicated by encircling with red lines. The geographic distribution of sublineage Ib from northernmost Honshu is shaded.

Fig. 2. Maximum likelihood tree of the mitochondrial cytochrome b gene (Cytb) sequences (1,140 bp) from Microtus montebelli (A). Node numbers indicate bootstrap values (1,000 replicates). Bar indicates number of substitutions per site. Bayesian phylogenetic tree with Cytb sequences (1,140 bp) from M. montebelli, showing means of divergence times (Mya) of nodes and horizontal bars of the interval of 95% HPDs estimated using BEAST (B). Median-joining networks constructed with Cytb sequences (1,140 bp) from M. montebelli (C). The six major lineages (lineages I–VI) identified in the phylogenetic analysis are indicated, as well as sublineages (haplotype groups) Ia and Ib with star-like patterns. The numbers of mutations between haplotypes are indicated by hatching. The sizes of circles are proportional to the sample size.

We estimated the divergence times among the six distinct lineages using BEAST with an evolutionary rate of 0.029 substitutions/site/million years (Myr) and the TN93 + G + I nucleotide substitution model. The results showed that the mean values and highest posterior density intervals (HPDs) for the ages of the five nodes representing the occurrence of the six lineages ranged from 161,000 years ago (95% HPD, 9,300–241,000 years ago) to 308,000 years ago (95% HPD, 205,000–437,000 years ago) (Fig. 2B).

We constructed a median-joining (MJ) network using the Cytb sequences of M. montebelli (Fig. 2C). Lineage I formed a star-shaped-like configuration in its overall structure. In addition, it encompassed two star-shaped nested clusters (haplotype groups), designated as Ia and Ib. The “Sado cluster” (Ia) comprised eight haplotypes from Sado Island and two from northern Honshu. The haplotype group Ib consisted primarily of haplotypes from northernmost Honshu. Tajima’s D and Fu’s F_S mostly supported the hypothesis of sudden population expansion in lineage I and haplotype groups Ia and Ib (Table 1).

Table 1. Expansion time estimation for haplotype groups of Microtus montebelli and Lepus brachyurus

Species	Haplotype group	Main geographic areas	n	Mean pairwise differences (SD)	Tajima's D	Fu's FS	SSD (P)	RI (P)	t^a (5%–95%)	Estimated time with three evolutionary rates^b
Species	Haplotype group	Main geographic areas	n	Mean pairwise differences (SD)	Tajima's D	Fu's FS	SSD (P)	RI (P)	t^a (5%–95%)	0.028	0.047	0.11
Microtus montebelli
	I	N Honshu	33	3.53	-1.94**	-22.15**	0.001	0.024	3.60	56,300	33,600	14,400
				(1.84)			(0.76)	(0.60)	(2.52–4.56)	(39,500–71,400)	(23,500–42,600)	(10,000–18,200)
	Ia	Sado	10	2.02	-1.22	-5.03**	0.052	0.179	2.32	36,300	21,700	9,300
				(1.23)			(0.12)	(0.15)	(1.27–3.74)	(19,900–58,600)	(11,900–34,900)	(5,100–14,900)
	Ib	Northernmost Honshu	9	0.89	-1.61*	-2.81**	0.038	0.224	1.13	17,700	10,500	4,500
				(0.68)			(0.26)	(0.21)	(0.00–1.97)	(0–30,900)	(0–18,400)	(0–7,900)
Lepus brachyurus
	N2	N Honshu	47	3.84	-2.05**	-25.95**	0.005	0.032	3.93	61,600	36,700	15,700
				(1.99)			(0.22)	(0.17)	(2.85–4.56)	(44,600–71,400)	(26,600–42,500)	(11,300–18,100)
	N2a	C Honshu	12	1.44	-1.48	-3.66**	0.006	0.097	1.65	25,800	15,400	6,600
				(0.94)			(0.69)	(0.53)	(0.74–2.74)	(1,200–42,900)	(6,900–25,600)	(3,000–10,900)
	N2b	Sado	11	2.25	-1.12	-4.01**	0.045	0.150	2.77	43,400	25,800	11,000
				(1.34)			(0.50)	(0.13)	(1.57–3.96)	(24,600–62,000)	(14,600–37,000)	(6,300–15,800)
	N5	C & W Honshu	75	4.03	-1.99**	-25.9**	0.002	0.020	4.31	67,500	40,200	17,200
				(2.03)			(0.35)	(0.44)	(2.67–5.16)	(41,800–80,800)	(25,500–48,200)	(10,600–20,600)
	N5a	W Honshu	13	2.08	-1.67*	-3.5*	0.006	0.032	2.77	43,400	25,800	11,000
				(1.24)			(0.77)	(0.95)	(0.54–4.70)	(8,500–73,600)	(5,000–43,900)	(2,200–18,700)
	S1a	W Honshu	23	1.16	-1.73*	-3.95**	0.003	0.037	1.48	23,200	13,800	5,900
				(0.78)			(0.78)	(0.92)	(0.22–2.60)	(3,300–40,700)	(2,100–24,300)	(900–10,400)

n, number of sequences analyzed; SD, standard deviation; SSD, sum of squared deviations; RI, raggedness index ; *, P < 0.05; **, P < 0.02. Values of SSD, RI, and t under the demographic expansion model in the Arlequin analysis are presented.

^aLower and upper 95% confidence interval limits are shown in parentheses.

^bTime (years ago) when expansion event started was estimated using the t value obtained and three evolutionary rates (substitutions/site/million years). The shaded dates are based on the recommended rates proposed by Suzuki (2021).

The τ value for the entire sample set of lineage I was 3.60 (Table 1). Within this lineage, the internal haplotype groups, Ia and Ib, had τ values of 2.32 and 1.13, respectively. Using the predicted μ value of 0.11 substitutions/site/Myr, the onset of expansion for the entire lineage I was estimated to be 14,400 years ago (confidence interval (CI): 10,000–18,200 years ago). The haplotype groups Ia and Ib were estimated to have expanded 9,300 years ago (CI: 5,100–14,900 years ago) and 4,500 years ago (CI: 0–7,900 years ago), respectively.

Cytb sequence variation in L. brachyurus

The construction of an ML tree using the Cytb sequence data of L. brachyurus (n = 226, 1,140 bp) led to the identification of two distinct clades, N and S, mainly observed in northern and central Honshu and southern Honshu/Shikoku/Kyushu (Supplementary Fig. S1A, S1B), as described previously (Nunome et al., 2010). As demonstrated previously, lineages N and S consisted of five and two sublineages designated as N1–N5 and S1–S2, respectively. The monophyly of the haplotype groups N4 and S1 was supported with moderate bootstrap values of 84% and 69%, respectively.

From a BEAST analysis, the divergence between lineages N and S was estimated to be approximately 702,000 years ago, with a 95% HPD range of 504,000–891,000 years ago (Supplementary Fig. S1C). The times of the most recent common ancestor (tMRCA) of lineages N and S were estimated to be 208,000 years ago (95% HPD: 131,000–262,000 years ago) and 210,000 years ago (95% HPD: 132,000–269,000 years ago), respectively.

The constructed MJ network demonstrated that multiple sublineages within both the N and S lineages exhibited a star-like configuration, indicative of secondary population growth events (Fig. 3). The N2 sublineage was observed to comprise two internal clusters of a star-shaped pattern, designated as N2a and N2b. In consideration of the observed tendency of the N2b cluster to align with Sado-derived haplotypes, a Sado-derived haplotype (Fig. 3, asterisk) was incorporated into the group for comparative purposes, despite its potential status as an outlier. The N5 sublineage exhibited the emergence of 12 branches from a central ancestral haplotype, many of which subsequently underwent consolidation to form a haplotype group (N5a–N5h). One of these branches, N5a, which contains a relatively large number of haplotypes, exhibited a distinctive star-shaped configuration. While the majority of the N5 sublineage is constituted by haplotypes that originated from eastern Honshu, some groups were composed of haplotypes from western Honshu, the Kii Peninsula (N5a, N5f) and the Chugoku region (N5c, N5g). Lineage S exhibited grouping of the two sublineages S1 and S2. The haplotypes that constituted S1 were widely distributed throughout the Chugoku, Shikoku and Kyushu regions. The lineage exhibited an internal cluster of haplotypes (S1a) displaying a distinct star-shaped pattern, comprising haplotypes originating from the Chugoku region, western Honshu. The haplotypes that comprised S2 were primarily derived from Kyushu, Shikoku and the Oki Islands.

Fig. 3. A median-joining network constructed with mitochondrial cytochrome gene b sequences (1,140 bp) from Lepus brachyurus. The five sublineages of lineage N (N1–N5) and two sublineages of lineage S (S1, S2) are marked with shaded boxes. The haplotype groups exhibiting a star-shaped pattern in the N1, N2 and S1 sublineages (N2a, N2b, N5a and S1a) are indicated by bold dotted lines. The N2b haplotype group is composed of all Sado-derived haplotypes, including one located at a slightly distant position (asterisk). The remaining N5 sublineages are indicated by thin dotted lines, and the names of the regions are indicated. Mutation numbers between haplotypes are indicated by hatching.

We performed neutrality tests and mismatch distribution analysis on haplotype groups that showed a distinctive star-shaped pattern in the haplotype network of L. brachyurus. Six of the haplotype groups listed in Table 1 showed a rapid population expansion signal. Using an evolutionary rate of 0.11 substitutions/site/Myr, the onsets of the population expansion events of the two sublineages N2 and N5 and the two internal haplotype groups N2b and N5a were estimated to be ~15,000 years ago and ~10,000 years ago, respectively. The haplotype groups in central Honshu (N2a) and the Chugoku region (S1a) were estimated to have expanded relatively recently at 6,600 years ago (CI: 3,000–10,900 years ago) and 5,900 years ago (CI: 900–10,400 years ago), respectively.

DISCUSSION

Influence of global climatic changes on genetic structure and demographic history

We examined the mtDNA diversity of the vole species M. montebelli, whose genetic structure has not been investigated previously. Phylogenetic analysis provisionally identified six distinct lineages, four of which have a geographically partitioned distribution spanning from northern Honshu to Kyushu (Fig. 1). The branching patterns of these six lineages could not be fully resolved (Fig. 2A), implying rapid divergence over a short evolutionary time. The six lineages were estimated to have diverged between approximately 200,000 and 300,000 years ago, suggesting the involvement of paleoclimatic abrupt warming (~240,000 years ago) following the antepenultimate glacial maximum (APGM) during MIS 8.

The two major haplotype groups of L. brachyurus, lineages N and S, were estimated to have diverged 702,000 years ago (95% HPD: 504,000–891,000 years ago) (Supplementary Fig. S1C). The results of BEAST analysis indicate that the sublineages of each of the lineages N and S diverged approximately 150,000 years ago. These two major lineages probably began to diversify during the period of rapid warming (~130,000 years ago) that followed the PGM in MIS 6. Based on the voles discussed above, it is reasonable to assume that the rapid expansion events triggered by the 100,000-year Quaternary glacial–interglacial cycle formed the basic framework for intraspecific variation in Japanese terrestrial mammals of temperate origin. Kai et al. (2024) reported rapid population expansion events in the black rat (Rattus rattus Complex) of East and Southeast Asia, both post-APGM and post-PGM. Therefore, traces of such rapid population expansion events in terrestrial small mammals inhabiting mid- and low-latitudes occurring at relatively old time points may tend to remain as intraspecific mtDNA diversity.

Two-stage rapid post-LGM expansion events

The MJ network had clusters with relatively small star-shaped patterns at the end of branches in lineages from northern and central Honshu in both species of voles and hares (Figs. 2C, 3). In voles, the τ value of lineage I was 3.60 (CI: 2.52–4.56), while the value for its internal haplotype group Ia was 2.32 (CI: 1.27–3.74) (Table 1). Similarly, the τ values of N2 and N5 of the hares showed relatively large τ values (3.93 and 4.31), with additional clusters N2b and N5a having smaller τ values (2.77). Notably, the presence of two large and small τ values in a single lineage provides strong evidence for two consecutive expansion events. The estimated expansion times of the large and small τ groups in both voles and hares were approximately 15,000 and 10,000 years ago, respectively. In consideration of paleoclimatic evidence, the two expansion events appear to correlate with the Bølling–Allerød and post-Younger Dryas abrupt warming periods (Schulz et al., 1998; Alley, 2000). These results are in accordance with the findings of previous studies on voles (genus Myodes), moles (genus Mogera) and wood mice (Honda et al., 2019; Nakamoto et al., 2021; Inoue et al., 2022).

The study provided two notable insights into the rapid expansion phenomena associated with the LGM. First, the effects of the LGM were found to exhibit a latitudinal gradient. The populations of the voles and hares were found to be significantly affected in northern and central Honshu after the LGM, whereas the effects were less pronounced in the southern and western regions. Second, the two-stage population expansion event after the LGM showed a distinctive spatial evolutionary pattern, with spatial expansion occurring in two stages. For example, in the N5 sublineages of hares, the ancestral haplotype and the secondary non-dispersed haplotype were distributed in northern parts (e.g., Hokuriku and Kanto regions) of the N5 range, indicating that N5 may have originated in the northern part of Honshu. The group of haplotypes that underwent the second stage of population expansion, which appears to have occurred in a novel space, were distributed not only in the central part of the distribution (e.g., Tokai region), but also in the Kansai and Chugoku regions (see Supplementary Fig. S1 and Fig. 3). The results indicate that the two expansions represent the spread of a single mtDNA lineage from the original area of the expansion event to the peripheral areas. This process, postulated by previous studies on Japanese small rodents (Tomozawa and Suzuki, 2008; Inoue et al., 2022), is of particular significance with regard to the geographic structuring of mtDNA.

In this study, we detected the presence of two haplotype groups with a star-shaped pattern with substantially low τ values that evolved in localized areas in voles and hares: the haplotype group Ib of M. montebelli (τ = 1.13) is distributed in northern Honshu (Fig. 2A) and the N2a and S2b sublineages of L. brachyurus (τ = 1.65) are distributed mainly in the Kanto and Chugoku regions, respectively (Supplementary Fig. S1A). Although the confidence interval based on an evolutionary rate of 0.11 substitutions/site/Myr is large and the possibility of a dispersal event after the Younger Dryas period cannot be excluded, it is possible that expansion events occurred 4,000 to 7,000 years ago (Table 1). Nakamoto et al. (2021) analyzed mtDNA of the large Japanese mole (Mogera wogura) and found evidence of population expansion in central Honshu about 5,000 to 6,000 years ago, implying that the expansion of plains after the end of the Jomon transgression was involved in the mole expansion event in the middle Holocene. Further studies are needed to increase the population size and address this issue in greater detail to determine whether these events were synchronized at some point during the middle Holocene.

Migration across the Sado Strait driven by the rapid expansion event

The results of the network analysis indicated that the haplotype groups, including those derived from the Sado Island population, exhibited a star-shaped pattern in both species of voles and hares examined in this study. The application of the elevated evolutionary rate (0.11 substitutions/site/Myr) yielded an estimated age of approximately 15,000 years ago for these expansion events in both species (Table 1). Haplotype clusters, predominantly those derived from Sado Island, were identified as internal clusters in both species. Their estimated age of occurrence was approximately 10,000 years ago. It seems reasonable to assume that the mtDNA haplotype of the Sado Island ancestor was introduced to Sado from Honshu at a minimum of 10,000 years ago and a maximum of 15,000 years ago. It seems reasonable to posit that this introduction was associated with the rapid population expansion events that occurred after the LGM. During the LGM, the northern limit of beech, which represents temperate broad-leaved trees, was located at 38°N in the Japan Sea coastal lowlands (Momohara and Ito, 2023). It is conceivable that this area served as a refuge for temperate mammal species. During the post-glacial warm period, temperate broad-leaved trees developed, and their populations may have expanded rapidly. In response to these changes, the population size of animals on Honshu may have increased rapidly during the warm period, and the habitat for animals may also have improved on Sado Island. We hypothesize that these conditions drove the migration to Sado. Furthermore, the arrival of the Bølling–Allerød interglacial and the warm period of the Holocene may have facilitated the transfer to Sado, as the land masses of Honshu and Sado were closer together than they are today (Bintanja et al., 2005). Further verification by analysis of other terrestrial mammals on Sado Island, such as moles and shrews, will be necessary to confirm this hypothesis.

CONCLUSION

We compared mtDNA variation in two temperate herbivore species, M. montebelli and L. brachyurus. Both species had multiple regional population lineages that tended to diverge simultaneously. The divergence of the two major regional lineages of voles and hares is assumed to have originated from expansion events following past glacial maxima, post-APGM and post-PGM, respectively. The effect of the LGM varies with latitude, with stronger effects in northern and weaker in southern populations. Following the LGM, population expansion was characterized by two successive phases, estimated to have occurred about 15,000 and 10,000 years ago. It has been demonstrated that these two consecutive rapid population expansions facilitated the geographic expansion of descendant lineages of specific haplotype groups. Furthermore, the expansion of this post-LGM population seems to have played a role in the dispersal and settlement of mtDNA in the Sado Strait in the vole and hare species. It would be beneficial for future research on the voles and hares to employ various mtDNA genes and to investigate genome-wide nucleotide variation, which have been demonstrated to facilitate a more comprehensive understanding of geographic structure, population dynamics and gene flow of the species (e.g., Fujiwara et al., 2024; Kinoshita et al., 2024; Ohdachi et al., 2024).

MATERIALS AND METHODS

Samples

Fifty samples of M. montebelli were collected from Honshu and Kyushu. The study used tissue and fecal samples from 27 individuals of L. brachyurus collected from Sado Island, the Oki Islands and the Chugoku region.

DNA extraction and sequence analysis

For M. montebelli, DNA was extracted using a DNeasy Blood & Tissue Kit (QIAGEN, Hilden, Germany). All DNA extracts were deposited at the Botanic Garden & Museum, Field Science Center for Northern Biosphere, Hokkaido University, Sapporo (HUNHM; Supplementary Table S1). Amplification of Cytb (1,140 bp) was performed by polymerase chain reaction (PCR). Primers were designed to target the 5' half Cytb fragment with L14115 (5'-GACATGAAAAATCATCGTTG-3'; Yasuda et al., 2005) and Mmon-R1 (5'-GAAGTCTTTGATTGTATAGTAGG-3'), and the 3' half with Mmon-F2 (5'-ACCCTAGTAGAATGAATCTGAGG-3') and H15300 (5'-GTTTACAAGACCAGAGTAAT-3'; Yasuda et al., 2005). PCR was performed with an initial denaturation step of 10 min at 95 °C followed by 30 cycles of 30 s at 95 °C, 30 s at 50 °C and 60 s at 60 °C, with a final extension step for 7 min at 72 °C, using an Ampli Taq Gold 360 Master Mix kit (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA). The PCR products were sequenced using a BigDye Terminator Cycle Sequencing kit v3.1 (Applied Biosystems, Thermo Fisher Scientific, Waltham, MA, USA) and an ABI3130 automated sequencer (Applied Biosystems). The sequences were edited with MEGA X (Kumar et al., 2018). For hares, DNA extraction and determination of the Cytb sequences were performed as described previously (Nunome et al., 2010).

Phylogenetic analysis and estimation of divergence times

Phylogenetic trees were constructed based on mtDNA sequences using the ML method implemented in MEGA X (Kumar et al., 2018). The Cytb data set (1,140 bp) for M. montebelli consisted of the 50 newly determined sequences and three sequences from databases (MF099543, 1,093 bp; NC_059928, 1,140 bp; AF163900, 1,140 bp), with Microtus oeconomus (AY220019) included as an outgroup (Supplementary Table S1). The L. brachyurus data set consisted of 199 sequences examined previously (Nunome et al., 2010), 27 sequences determined in this study (Supplementary Table S2) and two sequences of the outgroup taxon Lepus timidus (AJ279424 and AJ279425). A nucleotide substitution model was selected based on the Akaike Information Criterion using MEGA X (Kumar et al., 2018). The TN93 + G + I model was identified as the best and second-best fit model for the data sets of M. montebelli and L. brachyurus, respectively. The nodal support was assessed by 1,000 bootstrap replicates to estimate the reliability of the inferred phylogenetic relationships.

A Bayesian analysis of the branching times for both species was conducted using BEAST version 1.8.4 (Drummond and Rambaut, 2007), with the TN93 + G + I substitution model, tree prior distribution set to “coalescent; constant size,” and following previously described methods (Honda et al., 2019). We estimated the branching points of the species to have occurred more than 100,000 years ago. The evolutionary rates of p-distance in rodents, evaluated by Suzuki (2021), were proposed to be approximately 0.028 substitutions per site per Myr. This rate has been reported to be roughly equivalent in rodents and moles (Nakamoto et al., 2021), and was employed in this study for the analysis of the voles and hares. A comparison of the base substitutions of the TN93 + G substitution model with those of p-distance, as available in MEGA X, revealed that the former were approximately 104% and 108% of the latter, respectively. Consequently, evolutionary rates of 0.029 and 0.030 substitutions per site per year were applied for M. montebelli and L. brachyurus, respectively, for the BEAST analysis using the TN93 + G + I substitution model.

Demographic analysis

We constructed MJ haplotype networks using PopArt (Leigh and Bryant, 2015) with the default settings. We focused on star-shaped clusters that imply a sudden population expansion. For each of these clusters, we conducted Tajima’s D and Fu’s F_S tests of neutrality and calculated nucleotide diversity using Arlequin ver. 3.5.1.2, with significance levels set to P < 0.05 and P < 0.02, respectively (Excoffier and Lischer, 2010). We performed mismatch distribution analysis to examine the pattern of unimodality, which is expected in the case of a sudden expansion. The expected distribution was simulated under the sudden expansion model (Schneider and Excoffier, 1999; Excoffier, 2004) and evaluated using a parametric bootstrap approach (1,000 replicates). To examine deviations from the model of population expansion, we used the sum of squares of the deviation (SSD) goodness-of-fit test in Arlequin ver. 3.5.1.2. As a test statistic for the estimated sudden expansion models, we used the raggedness index (r; Harpending, 1994). Non-significant SSD and r values (P > 0.05) indicate the occurrence of a demographic expansion event. Mismatch distribution analysis provides τ as a statistical value proportional to the time since expansion. The CI around the estimated τ parameter was obtained through 1,000 replicates of coalescent simulation. The time in generations (t) since expansion of each cluster was estimated using the formula t = τ/2u, where u is the substitution rate per generation per haplotype sequence (Rogers and Harpending, 1992). We used the formula T = tg to calculate the time in years (T) since the expansion, where g represents the number of generations per year. In this context, the evolutionary rate (u) can be expressed as u = μgk, where μ represents the rate of substitutions per site per year and k denotes the sequence length of the haplotype. Applying these considerations (T = τg/2μgk), we derived the simplified equation T = τ/2μk, in which the information pertaining to the generation time has been canceled. The molecular clock rate for intraspecific differentiation has been shown to be time-dependent and to decrease with increasing time scale (Ho et al., 2011). Based on the findings of previous studies (Suzuki et al., 2015; Hanazaki et al., 2017; Nakamoto et al., 2021), we categorized the obtained τ values into three groups: < 4, 5–6 and > 7. For each category, we used specific evolutionary rates of 0.11, 0.047 and 0.028 substitutions/site/Myr, respectively.

COMPETING INTERESTS

The authors have no competing interests to declare.

AUTHOR CONTRIBUTIONS

H. S., T. Y., M. N. and G. K. conceived and designed the study. H. S., T. Y., T. E., M. H. and G. K. performed the field study. H. S., T. Y. and G. K. determined sequences and analyzed the data. H. S., N. M. and G. K. wrote the original draft of the manuscript. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS

The authors are grateful to the late Kimiyuki Tsuchiya for his collaborative research efforts in this study. We thank Masaru Kato, Satoshi Ohdachi and Shumpei Yasuda for their support in conducting this study. This study was supported by JSPS KAKENHI Grant Number JP18H05508 and partly supported by the Grant for Joint Research Program of the Institute of Low Temperature Science, Hokkaido University (24G031).

REFERENCES

Corresponding author

Register with J-STAGE for free!