- Research Article
- Open Access
Estimating the molecular evolutionary rates of mitochondrial genes referring to Quaternary ice age events with inferred population expansions and dispersals in Japanese Apodemus
BMC Evolutionary Biology volume 15, Article number: 187 (2015)
Determining reliable evolutionary rates of molecular markers is essential in illustrating historical episodes with phylogenetic inferences. Although emerging evidence has suggested a high evolutionary rate for intraspecific genetic variation, it is unclear how long such high evolutionary rates persist because a recent calibration point is rarely available. Other than using fossil evidence, it is possible to estimate evolutionary rates by relying on the well-established temporal framework of the Quaternary glacial cycles that would likely have promoted both rapid expansion events and interisland dispersal events.
We examined mitochondrial cytochrome b (Cytb) and control region (CR) gene sequences in two Japanese wood mouse species, Apodemus argenteus and A. speciosus, of temperate origin and found signs of rapid expansion in the population from Hokkaido, the northern island of Japan. Assuming that global warming after the last glacial period 7–10 thousand years before present (kyr BP) was associated with the expansion, the evolutionary rates (sites per million years, myr) of Cytb and CR were estimated as 11–16 % and 22–32 %, respectively, for A. argenteus, and 12–17 % and 17–24 %, respectively, for A. speciosus. Additionally, the significant signature of rapid expansion detected in the mtDNA sequences of A. speciosus from the remaining southern main islands, Honshu, Shikoku, and Kyushu, provided an estimated Cytb evolutionary rate of 3.1 %/site/myr under the assumption of a postglacial population expansion event long ago, most probably at 130 kyr BP. Bayesian analyses using the higher evolutionary rate of 11–17 %/site/myr for Cytb supported the recent demographic or divergence events associated with the Last Glacial Maximum. However, the slower evolutionary rate of 3.1 %/site/myr would be reasonable for several divergence events that were associated with glacial periods older than 130 kyr BP.
The faster and slower evolutionary rates of Cytb can account for divergences associated with the last and earlier glacial maxima, respectively, in the phylogenetic inference of murine rodents. The elevated evolutionary rate seemed to decline within 100,000 years.
To investigate the recent evolutionary history of organisms, many phylogeographic studies have assessed molecular dating using rapidly evolving markers, such as mitochondrial DNA (mtDNA) gene sequences. However, a gap has recently been suggested to exist in the evolutionary rates of mtDNA gene sequences between interspecific divergence and intraspecific polymorphism and, therefore, the assumption that a single molecular clock applies throughout all lineages over time may not be correct (e.g. [1–4]). In fact, to analyze evolutionary history on a recent timescale, substantially higher rates [e.g., 38.9 or 40 %/site/million years (myr)] in the mtDNA control region (CR) have been proposed in studies of rodents such as the field vole (Microtus agrestis ) and house mouse (Mus musculus ), based on paleoclimatological and/or geological events using Bayesian inferences, contrary to the rates “traditionally” used for mammals (e.g., 2–4 %/site/myr; [7, 8]). However, how long such a high evolutionary rate is applicable [9, 10] remains unclear due to the lack of reliable calibration points, such as fossil evidence below the species level.
Apodemus species in the Japanese Islands are ideal subjects for addressing this issue. In Eurasia, wood mice (Apodemus spp.) are widely distributed throughout temperate and subarctic regions, and their evolution has been spurred by habitat topography and the expansion of the temperate zone in the Tertiary period [11–14]. Two endemic Apodemus species in Japan (A. argenteus and A. speciosus) are known to have a long (5–6 myr) evolutionary history, independent from that of congeneric continental species according to molecular phylogenetic studies [14, 15]. Because both of these species are common in temperate habitats, such as broadleaf forests, their demography appears to be associated with the reduction and expansion of temperate environments during the Quaternary glacial cycles, as suggested in previous phylogeographic studies [16, 17]. Thus, one can estimate a molecular evolutionary rate by associating a specific glacial cycle with a sign of population size change in a phylogenetic diversification pattern. Recent studies on pollen stratigraphy, for example, revealed that the reduction in temperate broadleaf forest during the Last Glacial Maximum (LGM) was severe in the northern Japanese Islands, such as Hokkaido and northern Honshu [18, 19]. Thus, we may be able to use this event as a calibration point to estimate a molecular evolutionary rate over a short timescale.
Another advantage of studying these two Apodemus species is that we can use information on the timing of gene flow among isolated island populations to evaluate the estimated evolutionary rates. Both species inhabit almost all Japanese Islands. The Japanese Islands comprise the four main islands of Hokkaido, Honshu, Shikoku, and Kyushu arranged from northeast to southwest, as well as several insular regions, such as Sado Island, the Izu Islands, the Tsushima Islands, and the Satsunan Islands. Some of these islands are separated from the main islands by deep sea straits (depth 100–200 m) that seem to have limited the dispersal of terrestrial animals except when the sea level dropped during Quaternary glacial periods. In particular, the populations of A. speciosus on peripheral islands will provide useful information. For example, the divergence time between the Hokkaido and Honshu populations of A. speciosus is similar to those of other island populations isolated by deep sea straits, suggesting that colonization events into the peripheral islands occurred simultaneously when the sea level dropped in the Quaternary glacial periods [16, 17, 20]. Such a priori information about periodic gene flow to isolated island populations will help in evaluating the validity of the evolutionary rate estimates.
In this study, we focused on the postglacial population expansion of the two Apodemus species in Hokkaido, and then estimated the evolutionary rates of two mtDNA sequences, cytochrome b (Cytb) and CR, the most popular markers for phylogenetic and phylogeographic inferences in mammals. One can reasonably presume that the two species, which are of temperate origin, would have been affected by the Quaternary glacial periods , and that the population in Hokkaido, the northernmost main island of Japan, would have been especially impacted. We also evaluated the validity of the estimated evolutionary rates by comparing the sequence divergence time calculated using the evolutionary rates and the temporal aspects of the divergence events between the four main islands and peripheral islands, which are associated with periodic changes in sea level in the straits.
In total, we used 134 A. argenteus individuals (58 localities) and 128 A. speciosus individuals (51 localities) (Fig. 1, Table 1). The samples included DNA and tissue samples stored in our laboratory [14, 16], in addition to new specimens from Hokkaido, collected between 2011 and 2013. These animals were captured specifically for this study, and appropriate permissions were obtained from all relevant prefectural authorities. All experiments involving the sacrifice of live animals captured in the field were conducted within the ethical guidelines of the Mammal Society of Japan.
DNA was extracted using a standard phenol-chloroform extraction method . For both species, two mtDNA sequences were analyzed: the entire coding region of Cytb (1,140 bp) and half (the hypervariable region I) of CR (A. argenteus, 560 bp; A. speciosus, 554 bp). Sequence determination for Cytb was performed with semi-nested polymerase chain reaction (PCR), amplifying two fragments (“upper” and “lower”) in the second PCR . The first run was done using the primers L14724 and H15915 , at a final concentration of 0.05 pmol/μl, and Amplitaq-Gold DNA polymerase (Applied Biosystems) under thermal cycling parameters of 95 °C for 10 min and 30 cycles at 95 °C for 30 s., 50 °C for 30 s., and 60 °C for 30 s. The second PCR was done with primers at 0.05 pmol/μl for both upper (R-L14724: 5′-CAGGAAACAGCTATGACCGATATGAAAAACCATCGTTG-3′, SNH655: 5′-TGTAAAACGACGGCCAGTTGTGTAGTATGGGTGGAATGG-3′) and lower (SNL497: 5′-CAGGAAACAGCTATGACC CCTAGTAGAATGAATCTGAGG-3′, R-H15916: 5′-TGTAAAACGACGGCCAGTGTCATCTCCGGTTTACAAGA-3′) regions under conditions of 30 cycles at 96 °C for 30 s, 50 °C for 60 s, and 60 °C for 30 s, using Amplitaq DNA polymerase (Applied Biosystems). PCR for the CR sequences was performed with the AmpliTaq Gold 360 Master Mix kit (Applied Biosystems) using primers L15399 (5′-GCACCCAAAGCTGATATTCT-3′; ) and CR1 (5′-CATGCCTTGACGGCTATGTT-3′; ) at a final concentration of 0.5 pmol/μl, using PCR conditions of 96 °C for 10 min and 35 cycles at 96 °C for 30 s, 50 °C for 60 s, and 72 °C for 60 s. The PCR products were sequenced using the PRISM Ready Reaction DyeDeoxy Terminator Cycle Sequencing Kit (Applied Biosystems) and an ABI3130 automated sequencer (Applied Biosystems). The sequences of both strands were determined using the universal primers (M13RP1 and −21 M13; Applied Biosystems) in Cytb and the primers used for PCR in CR. The sequences were aligned using ClustalW implemented in MEGA5 .
Phylogenetic trees with the concatenated mtDNA sequences of Cytb and CR were generated with MEGA5 , using maximum-likelihood (ML), maximum-parsimony (MP), and neighbor-joining (NJ). All indels were excluded from the phylogenetic analysis. In addition to sequences of Apodemus peninsulae (HS123; collected from Hokkaido), A. speciosus (HS97) and A. argenteus (HS223) were used as outgroup taxa in the phylogenetic analyses of A. argenteus and A. speciosus, respectively. In the ML analysis of the A. argenteus data, the Tamura three-parameter model (TN92 + G + I) was selected using MEGA5 with likelihood-ratio tests and the Akaike information criterion (AIC). The Hasegawa–Kishino–Yano model (HKY85) was used in the ML tree of A. speciosus. In the NJ analysis, the Kimura two-parameter model (K2) was used. The reliability of nodes was assessed using 1000 bootstrap replicates. Median-joining networks were constructed for both markers using the Network program (ver. 18.104.22.168) . We used ARLEQUIN (ver. 3.5) to calculate the nucleotide (π) and haplotype (Hd) diversity using all data .
To detect rapid population expansion, mismatch distribution analyses, neutrality tests (Tajima’s D  and Fu’s F S ) were performed using ARLEQUIN (ver. 3.5). The significance of neutrality was tested with 1,000 replicates of coalescent simulation. The mismatch distributions of the Cytb sequence were calculated for all data and for each phylogroup of both species. The mismatch distributions of the CR sequences were calculated for all populations and the Hokkaido populations of each species. A smooth, unimodal distribution indicates population growth . The expected distribution was simulated under the sudden expansion model [32, 33]. We tested the validity of the sudden expansion model using a parametric bootstrap approach with 1,000 replicates. In this method, for each replicate, the sum of squared deviation (SSD) between the observed and expected distributions is compared with the SSD between the simulated and expected distributions using ARLEQUIN (ver. 3.5). We used the raggedness index (r)  as a test statistic for the estimated sudden expansion models.
The expansion parameter tau (τ) was estimated using ARLEQUIN (ver. 3.5)  in each cluster in which signs of sudden demographic expansion were evident. The evolutionary rates of the mtDNA markers were estimated using the formula t = τ/2u, where t is the time since expansion (in generations) and u is the cumulative evolutionary rate per generation for the whole sequence . The value of u was derived from the formula u = μkg, where μ is the evolutionary rate per site per year, k is the sequence length, and g is the generation time in years. We assumed the average generation times of both species were 1 year because both species are believed to breed once or twice per year. The time since expansion (t) for the Hokkaido populations of both species was assumed to be 7–10 kyr BP, based on the timing of the recovery of the Quercus population in northern Hokkaido from the impact of the last glacial period .
The time to the most recent common ancestor (TMRCA) of the Cytb sequences and 95 % highest posterior density (HPD) were estimated using BEAST (ver. 1.7.5; http://beast.bio.ed.ac.uk) . The analysis was carried out with the HKY85 substitution model and the strict-clock model using expansion growth as the tree prior. Two European species of Apodemus sylvaticus (accession no. JF819979) and Apodemus flavicollis (accession no. JF819969) whose divergence time was estimated at 2–3 million years ago , were used as outgroup taxa. Bayesian searches were conducted using the Markov chain Monte Carlo (MCMC) method for 10 million generations. The first 1 million generations were discarded as a burn-in. Independence among samplings was confirmed for each run by checking for high effective sample sizes (>200). Tracer (ver. 1.5)  was used to assess convergence of the MCMC chains.
Sequence diversity indices for the two gene data sets are summarized in Table 2. Complete Cytb sequences (1,140 bp) were obtained from 134 and 128 specimens of A. argenteus (187 variable sites that defined 101 haplotypes) and A. speciosus (185 variable sites that defined 83 haplotypes), respectively. No deletion or insertion was detected in the Cytb sequences of either species. CR sequences for A. argenteus (560 bp; 134 specimens) and A. speciosus (554 bp; 123 specimens) contained 100 and 96 variable sites that defined 91 and 78 haplotypes, respectively. In A. argenteus, 11 indels were detected, whereas only one insertion was detected in A. speciosus.
Population genetic analyses with mtDNA
To determine the phylogeographic structures of the two Apodemus species, phylogenetic trees and median-joining (MJ) network trees were constructed using Cytb and CR sequences, and their concatenated sequences. The phylogenetic analyses with the Cytb sequences and concatenated sequences of Cytb and CR in A. argenteus (Fig. 2a) revealed that the haplotypes fell into two clusters (I and II), which were supported by low (<54 % in ML and MP) and high (>89 % in NJ) bootstrap values depending the methods used, but they were unsupported (<50 %) with the CR sequences alone (data not shown). The genetic distance for the Cytb sequences between the two major clusters was calculated as 2.9 % on average. Cluster I contained a relatively large subcluster, representing haplotypes from almost all Japanese islands, from Hokkaido to Kyushu, with low bootstrap values (<65 %). The phylogenetic analyses further revealed a cluster representing all the Hokkaido haplotypes (Ia-1) within subclade Ia, with low to moderate supporting values (64–84 %). Notably, a haplotype from eastern Honshu (HS158 from Gunma Prefecture) clearly was closely associated with the Hokkaido cluster Ia-1, although its supporting values were low (<67 %; Fig. 2a). Cluster II was subdivided into two clusters, IIa and IIb, with relatively high supporting values. Cluster IIa contained two haplotypes from the southernmost tip of eastern Honshu (Mt. Amagi, Shizuoka Prefecture; locality 38 in Fig. 1) and Kyushu (Miyazaki, locality 56), and Cluster IIb included eight haplotypes from Shikoku and Kyushu. The MJ networks of the Cytb, CR, and concatenated sequences (Fig. 3) showed star-like structures for the haplotypes from Hokkaido (Ia-1), as predicted in the phylogenetic tree analysis. A portion of haplotypes from both western and eastern Honshu also showed a star-like structure (Ia-2) in the Cytb data set, although the structure was ambiguous in the networks of CR and concatenated data sets.
The 83 haplotypes of the concatenated Cytb and CR sequences from A. speciosus fell into two major clusters (I and II) that were always supported, but to varying degrees depending on the methods used (Figs. 2b and 4). Cluster I included haplotypes from the three main islands of Japan (Honshu, Shikoku, and Kyushu) and their closely associated islands, such as the Oki and Tsushima Islands. The haplotypes from the Tsushima Islands were clustered together in cluster I. Cluster II encompassed four haplotype groups representing four remote island regions, Hokkaido (IIa), Sado Island (IIb), the Izu Islands (IIc), and the Satsunan Islands (IId), each of which was supported by the three tree-building methods, although the supporting values were all relatively low. The genetic distance of Cytb between the two clusters was 2.4 %. The MJ networks of the Cytb (Fig. 4a) and concatenated Cytb and CR sequences (Fig. 4c) showed that the clusters corresponded to those in the phylogenetic analysis of Cytb. The MJ network tree of the CR sequences showed no clear relationship (Fig. 4b), as was the case for A. argenteus. The networks of the Cytb and CR sequences were consistent with those from the phylogenetic tree and showed more defined relationships among the clusters (Figs. 2b and 4c). The MJ network constructed from the concatenated sequences revealed two distinct groups in the haplotypes from Hokkaido (IIa), with different levels of sequence divergence: IIa-1 (low divergence; π = 0.244 %) and IIa excluding IIa-1 (IIaex_IIa-1; high divergence; π = 0.525 %). The IIa-1 group was distributed throughout Hokkaido, whereas the IIaex_IIa-1 group was restricted to a few localities, such as Rishiri Island, Kunashiri Island, and the eastern part of Hokkaido. The Sado Island population (IIb) could be divided into subclusters IIb-1 and IIb-2, the Izu Islands populations (IIc) into subclusters IIc-1 and IIc-2, and the Satsunan Islands (IId) populations into subclusters IId-1 (“Osumi”) and IId-2 (“Tokara”), each of which was further subdivided into two geographic areas.
Genetic diversity indices, including the nucleotide diversity (π) and haplotype diversity (Hd) values, were calculated for each of the phylogroups denoted in this study (Table 2). For A. argenteus, the π values of the total samples were 1.091 % and 1.771 %, and the Hd values were 0.980 and 0.963, for Cytb and CR, respectively. Notably, the π values of both markers for the Hokkaido population were low (0.212 % and 0.416 %, respectively). For A. speciosus, the respective π values of the whole samples were 1.606 % and 1.426 %, and the Hd values were 0.990 and 0.991 for Cytb and CR, respectively. As with A. argenteus, the Hokkaido population showed quite low nucleotide and haplotype diversities. The Tajima’s D and Fu’s F S values for the Hokkaido population of A. argenteus were significantly negative for both markers. In A. speciosus, Tajima’s D and Fu’s F S values for the total sample set and the IIa-1 group were significantly negative. The IIaex_IIa-1 had significant negative values for Fu’s F S for both markers, although Tajima’s D was not significant. Some other clusters, such as the Ia-2 of A. argenteus (mainly those from the eastern part of Honshu) and cluster I of A. speciosus, also had significant negative values for the neutrality tests of both markers (Table 2).
The mismatch distribution analysis of the Cytb sequence was performed on each phylogroup except for two groups with small sample sizes (Sado Island and Izu Islands in A. speciosus): clusters I, II, Ia, Ia-1, and Ia-2 for A. argenteus, and clusters I, II, IIa, IId, IIa-1, and IIaex_IIa-1 for A. speciosus (Fig. 5). The mismatch distribution analysis with CR sequences was conducted only for the total samples and the Hokkaido population that showed clear clustering patterns. For the majority of the clusters studied, the sudden expansion model was not rejected with either marker, but was rejected for IIaex_IIa-1 of Cytb of A. speciosus (Table 3). In A. argenteus, analyses of the Hokkaido population showed unimodal distributions with similar τ values for both markers (τ = 2.547 and 2.482 for Cytb and CR, respectively). The Hokkaido population of A. speciosus tended to show bimodal distributions for both markers, although when separated, the two groups of Hokkaido haplotypes, IIa-1 and IIaex_IIa-1, had unimodal distributions for both markers, except Cytb from IIaex_IIa-1.
Evolutionary rate estimation
The evolutionary rates (μ, per site per myr) for both markers were calculated using the formula, t (time after expansion, in years) = τ/2μk (k: sequence length), suggesting that the Hokkaido populations experienced rapid postglacial expansion, with the assumption that the latest expansion occurred 7–10 kyr BP, in association with the expansion of the Quercus population in northern Hokkaido . We also accounted for the possibility that the expansion of the wood mouse population occurred immediately after the LGM (15 kyr BP), given that Apodemus includes habitat generalist species and that tree species in Hokkaido largely recovered after the end of the glacial maximum . The evolutionary rates of Cytb were estimated as 7.5–16.0 %/site/myr and 8.0–17.1 %/site/myr for the Hokkaido populations of A. argenteus (Ia-1) and A. speciosus (IIa-1), respectively (Table 3). The evolutionary rates of the CR sequence were 14.8–31.7 %/site/myr and 11.3–24.2 %/site/myr for A. argenteus and A. speciosus, respectively. We estimated the evolutionary rates of Cytb for the predicted expansion event of cluster I of A. speciosus occurring in the three main islands of Honshu, Shikoku, and Kyushu, with the assumption that the expansion started from the end of the penultimate glacial time, 130 kyr BP, and obtained a value of 3.1 %/site/myr (Table 3).
Divergence time estimation
The divergence times were estimated based on Bayesian inference using the evolutionary rates of Cytb obtained in the previous section. BEAST analysis was not performed with the CR sequences due to the greater inequality in branch lengths observed on the CR NJ trees compared with the Cytb NJ trees (data not shown), which suggested less regular substitution fixation rates over evolutionary time (e.g., ). In the BEAST analysis for A. argenteus (data not shown), the divergence time between the Hokkaido population and the closest related haplotype from Honshu (HS158/Gunma Prefecture) was estimated at 39.0, 34.6, and 22.2 kyr BP, using evolutionary rates of 7.5, 11.2 and 16.0 %/site/myr, respectively.
The other older divergences addressed in this study (node A for A. argenteus and nodes A’–H’ for A. speciosus; Figs. 2 and 6) were assessed using faster evolutionary rates of 6.7, 12.0, and 17.1 %/site/myr and the rather conservative evolutionary rates reported previously, 2.4 %/site/myr  and 2.8 %/site/myr , in addition to the value of 3.1 %/site/myr. The resulting estimated divergence times, including those of the two outgroup taxa, A. flavicollis and A. sylvaticus, are summarized in Table 4.
The mtDNA phylogeographic structures of Japanese wood mice of temperate origin showed signs of rapid population expansion and lineage diversification, presumably reflecting global climatic and sea-level changes during the Quaternary glacial cycles (Fig. 6). This would thus provide calibration points for assessing evolutionary rates and an opportunity for evaluating the validity of the predicted evolutionary rates and those obtained from interspecies phylogenetic inferences using calibration points from fossils.
Rapid expansion events of wood mice and associated climatic changes
Mismatch distribution analyses (Fig. 5) and neutrality tests (Table 2) suggested that both A. argenteus (τ = 2.55 for Cytb) and A. speciosus (τ = 2.73) expanded recently in Hokkaido. Considering the paleoenvironment of Hokkaido and the ecological features of the Apodemus species, the recent population expansions in Hokkaido most likely occurred coincidentally with rapid global warming after the LGM (~20 kyr BP). In northern Europe and North America, the present phylogeographic structures of temperate species have been influenced greatly by population size changes due to glacier extension in the LGM and the subsequent period of global warming . Similarly, the phylogeographic structures of temperate species in Japan were reportedly influenced by the LGM, especially in the northern part of the islands, including northern Japan [38, 39]. The impact of the LGM and subsequent Younger Dryas glacial re-advancement (YD, 11–12 kyr BP ) changed the tree species composition of Hokkaido forests from cold coniferous species to temperate broad-leaved species, including oak trees (Quercus) . Because both of the Apodemus species generally inhabit broad-leaved forests, the signatures of the recent expansion events of the two Apodemus species in Hokkaido are likely to have been associated with the expansion of the temperate broadleaf forests after the YD (~7–10 kyr BP ), although we can not exclude the idea that the expansion event occurred shortly after the LGM, 15 kyr BP.
The mismatch analysis and network analysis of Hokkaido A. speciosus populations (Figs. 4 and 5) suggested another expansion of IIaex_IIa-1, with a τ of 8.57. We observed a star-like pattern for this haplotype group in the networks with the CR and concatenate sequences but not in that for Ctyb. Given that the expansion event for IIa-1with τ = 2.73 is most likely to have occurred 7–15 kyr BP, the other expansion event for IIaex_IIa-1 is likely attributable to the next nearest global climate fluctuation (the rapid cooling and subsequent rapid warming): more specifically, the periods 60 kyr BP or 130 kyr BP, equivalent to the terminations of marine isotope stages 4 (MIS 4) and 6 (MIS 6), respectively (Fig. 6; ). Deciduous broad-leaved forests are known to have decreased and increased during these periods in Japan [18, 43]. We favor the latter, because the former glacial period was less severe, leading to an extreme bottleneck effect and because this view is consistent with the expansion event of cluster I of A. speciosus, which is most likely due to rapid warming at the end of MIS 6, as discussed later. The haplotypes belonging to the IIaex_IIa-1 group tended to show limited distributions on insular and peripheral parts of mainland Hokkaido near the coast, such as Rishiri Island, Kunashiri Island, and the eastern part of Hokkaido (Fig. 2b), indicating that in Hokkaido, the environment of the inland mountainous areas severely affected temperate species during the LGM, whereas coastal areas served as refugia. In contrast, A. argenteus was likely severely affected by the LGM, partly because of its preference for mountainous habitats.
The neutrality tests (Table 2) and mismatch distribution analysis (Fig. 5) suggested a historical population expansion of A. speciosus in Honshu, Shikoku, Kyushu, and some peripheral islands, such as the Tsushima Islands (cluster I, τ = 9.19 in Cytb). Comparing the τ value with the Hokkaido population (τ = 2.73), which is likely attributable to the LGM as discussed above, caused us to consider that the expansion of cluster I was likely associated with a glacial cycle older than the LGM and should have happened simultaneously or slightly before the earlier expansion presumed in Hokkaido (IIaex_IIa-1, τ = 8.57). Pronounced changes in regional vegetation at the beginning of the last interglacial period (128,000 years ago; MIS 5e) have been recorded in pollen-based quantitative biome reconstructions from the North Eurasian study sites  and central Honshu [18, 45]. Given that the postglacial appearance of Quercus trees is evident [18, 44, 46], the rapid expansion of A. speciosus (cluster I) in Honshu, Shikoku, and Kyushu is most likely due to climatic shifts in the last interglacial period. Notably, another rodent of broad-leaved forests, the Japanese giant flying squirrel Petaurista leucogenys, displays a similar trend of rapid expansion and mitochondrial nucleotide diversity in the northern part of its range; this expansion likely occurred at the onset of the last interglacial period . In addition, cluster I expansion generated the lineage of the Tsushima Islands (Fig. 4c); this lineage is separated from Kyushu by the Tsushima Strait, which is thought to have narrowed during the penultimate glacial maximum . Thus, the expansion of cluster I conceivably started 130 kyr BP, contemporaneous with the end of the penultimate glaciation.
Similarly, the mismatch distribution and neutrality test suggest glacial cycle-coupled expansion in cluster I, representing the main islands of Honshu, Shikoku, and Kyushu in A. argenteus (Fig. 5a, Table 2). By contrast, except for the Hokkaido subcluster, Ia-1, the network pattern of cluster I appears to have multiple subclusters (Fig. 3), indicating the involvement of multiple expansion events that are hard to assign to specific glacial cycles due to insufficient sample size. Cluster II consists of haplotypes from the southwestern part of the Japanese islands—Kyushu, Shikoku, and the southernmost tip of the Izu Peninsula (locality 38 in Fig. 1a)—suggesting effective refugia for anciently divergent mtDNA lineages during glacial periods. The presence of refugia along Honshu, Shikoku, and Kyushu has been suggested by phylogeographic studies on other mammals, such as hares , black bears , sika deer , flying squirrels , macaques , and voles .
A temporal view of the divergences among intraspecific lineages
Relying on the evolutionary history of the wood mice discussed above, we can estimate the molecular evolutionary rates of mtDNA over different timescales. First, attributing the latest population expansion in Hokkaido to the end of the last glacial period, the evolutionary rates of Cytb are estimated as 7.5–16 %/site/myr and 8.0–17.1 %/site/myr for A. argenteus and A. speciosus, respectively (Table 3). Second, if we accept that the earlier expansion of the Hokkaido population of A. speciosus started 130 kyr BP, the rate of Cytb is 2.9 %/site/myr. Finally, if we assume the expansion of A. specious in the southern islands of Honshu, Shikoku, and Kyushu started at 130 kyr BP, the evolutionary rate is 3.1 %/site/myr, which is comparable to those obtained previously (e.g., 2.8 % and 2.4 %) in the interspecies phylogenetic inferences with fossil evidence. Accordingly, these rate estimates can be categorized into faster (8–17 %) and slower (2.4–3.1 %) rates.
In A. speciosus, when we ran the BEAST analysis with the faster Cytb evolutionary rate, the resulting estimates for divergences were quite recent (e.g., 25–46 kyr BP for nodes C’–H’ with 17.1 %/site/myr; Table 4). For example, the divergences among the islands of Izu, Sado, and Hokkaido (node C’) and between Kyushu and the Satsunan Islands (node D’) were 46 and 41 kyr BP, respectively, suggesting dispersals in the MIS 3 period (30–57 kyr BP). However, this would seem highly unlikely because the sea-level drop at MIS 3 is thought to have been at most around −90 m (Fig. 6). Thus, most straits, including the Osumi Strait, which separates the Satsunan Islands and Kyushu, and the Tsugaru Strait between Honshu and Hokkaido, are supposed to have been too deep and wide in MIS 3 for interisland dispersals.
In contrast, the BEAST analysis with the lower evolutionary rate of 3.1 %/site/myr (Table 4, Fig. 6) suggests that the historical dispersals of A. speciosus to the remote islands happened at or immediately after either MIS 6 glacial periods (130–191 kyr BP) or older ones, of 100 kyr each [21, 53]. TMRCAs for two lineages of Izu (IIc-1, IIc-2; genetic distance), and Sado (IIb-1, IIb-2), Osumi (IId-1: Tanegashima vs. Yakushima and Kuchinoerabujima), and Tokara (IId-2: Kuchinoshima vs. Nakanoshima) islands were estimated to be 125 (kyr BP; 95 % HPD = 68–192.6 kyr BP), 131.9 (69.9–201.2), 157.5 (105–216.9), and 150.5 (86.4–219.5), respectively. Notably, the case of Sado Island is consistent with the geological evidence: at 110–130 kyr BP in the last interglacial period, Sado Island was separated into two islands by a sea strait . Additionally, the resulting estimates of TMRCAs of cluster I/II, cluster II, the lineages on Hokkaido/Sado/Izu, and Satsunan Islands are 463 (kyr BP; 95 % HPD = 365–569), 355 (252–478), 260 (190–333), and 229 (153–317), respectively (Table 4, Fig. 6). These timings may have been associated with the glacial periods of MIS 12 (424–478 kyr BP), MIS 10 (334–374 kyr BP), and MIS 8 (253–300 kyr BP), respectively. The sea-level drops during these periods were more than 100 m, and thus, the peripheral islands (Hokkaido, Sado, Satsunan, and Tsushima Islands) and the main islands were most likely connected or very close to each other.
Because the application of the faster evolutionary rate to relatively ancient divergences may lead to unrealistically underestimated values, as discussed above, we favor the TMRCA estimates of 589 kyr BP (95 % HPD = 453–729), 650 (500–809), and 750 (529–938) for clusters I and II of A. argenteus from the slower evolutionary rates of 3.1 %, 2.8 %, and 2.4 %, respectively (Table 4). Thus, it seems reasonable to propose that the ancient divergence was associated with the end of glacial periods, most likely that following MIS 16 (621–675 kyr BP). However, applying the faster evolutionary rate would be meaningful regarding a rather recent event: specifically, the divergence of the mtDNA between Hokkaido and Honshu in A. argenteus. From the BEAST analysis, the divergence time of the Hokkaido lineage from the most closely related haplotype of Honshu (Gunma, HS158; locality number 28 in Fig. 1a) was estimated as 22–39 kyr BP, with evolutionary rates of 7.5–16 %/site/myr. This indicates that the mtDNA of A. argenteus came to Hokkaido over the Tsugaru Strait (130 m) at most 39 kyr BP. The sea level of the Japan Sea during the LGM is presumed to have been ~120 m lower than that at present . If we take into account the trend of hydro-isostasy, namely that the depth of the Tsugaru Strait during the LGM was not as deep as currently seen due to absence of the weight of seawater (e.g., ), it then seems plausible to conclude that the recent migration of the mtDNA of A. argenteus occurred from Honshu to Hokkaido during or immediately after the LGM across a shallower, more narrow version of the Tsugaru Strait.
Time-dependent evolutionary rates for Apodemus species
Emerging evidence suggests that molecular evolutionary rates calculated over a few generations may differ substantially from those calculated by phylogenetic approaches [2, 9]. This is mainly because the molecular evolutionary rate measured over a short timescale counts mutations that are slightly deleterious and thus eventually eliminated from the population by natural selection (thus reflecting the mutation rate), whereas that calculated over a long timescale counts mutations that are fixed in the population by natural selection or genetic drift (thus reflecting the substitution rate) [4, 56]. However, it is unclear how long a high evolutionary rate remains over evolutionary timescales due to the lack of information about the time of sequence divergence at relatively long timescales, such as hundreds of thousands of generations .
As discussed, our data provide an ideal opportunity to assess the time dependency of molecular evolutionary rates over several hundred thousand years. The probable molecular evolutionary rates for recent timescales, such as 10 kyr BP, is faster; e.g., 8–17 %/site/myr and 11–32 %/site/myr for Cytb and CR, respectively (Table 3). In contrast, the use of the slower evolutionary rate (e.g., ~3 %/site/myr for Cytb) is more favorable for assessing older lineage divergences and population expansion events, 130 kyr BP or older. This in turn indicates that the evolutionary rate declines dramatically within 100 kyr. The presumed duration of the elevated evolutionary rate is shorter than those suggested for humans or insects (1–2 million years) [56, 57], and is comparable with that estimated in galaxiid fish (<200 kyr)  and cichlid fish (~100 kyr) . Future studies on other closely related species with different population histories may help to elucidate the time dependency of molecular rates and clarify the factors determining the duration of elevated evolutionary rates.
Overall, we show that the phylogeographic structures and population histories of two wood mouse species were greatly influenced by the Quaternary glacial cycles. Their population size changes and colonization events onto isolated islands were most likely associated with climatic and sea-level changes during the last and penultimate glacial cycles. Moreover, our work yielded fruitful results with respect to estimating the molecular evolutionary rates of mtDNA sequences (Cytb and CR), which change depending on the time that has passed over the last 100,000 years. These contributions are valuable in phylogeographic inferences, particularly for murine rodents.
Availability of supporting data
The nucleotide sequences reported in this paper appear in the DDBJ, EMBL, and GenBank nucleotide sequence databases under accession numbers AB974692–AB975183. Sequence data files in nexus file format, together with Supplementary Information files are available for downloading from the Dryad repository: doi:10.5061/dryad.4902q.
Ho SYW, Larson G. Molecular clocks: when timesare a-changin’. TRENDS Genet. 2006;22:79–83.
Howell N, Howell C, Elson JL. Time dependency of molecular rate estimates for mtDNA: this is not the time for wishful thinking. Heredity (Edinb). 2008;101:107–8.
Galtier N, Nabholz B, Glémin S, Hurst GDD. Mitochondrial DNA as a marker of molecular diversity: a reappraisal. Mol Ecol. 2009;18:4541–50.
Ho SYW, Lanfear R, Bromham L, Phillips MJ, Soubrier J, Rodrigo AG, et al. Time-dependent rates of molecular evolution. Mol Ecol. 2011;20:3087–101.
Herman JS, Searle JB. Post-glacial partitioning of mitochondrial genetic variation in the field vole. Proc R Soc B Biol Sci. 2011;278:3601–7.
Rajabi‐Maham H, Orth A, Bonhomme F. Phylogeography and postglacial expansion of Mus musculus domesticus inferred from mitochondrial DNA coalescent, from Iran to Europe. Mol Ecol. 2008;17:627–41.
Brown WM, George M, Wilson AC. Rapid evolution of animal mitochondrial DNA. Proc Natl Acad Sci. 1979;76:1967–71.
Wilson AC, Cann RL, Carr SM, George M, Gyllensten UB, Helm-Bychowski KM, et al. Mitochondrial DNA and two perspectives on evolutionary genetics. Biol J Linn Soc. 1985;26:375–400.
Ho SYW, Shapiro B, Phillips MJ, Cooper A, Drummond AJ. Evidence for time dependency of molecular rate estimates. Syst Biol. 2007;56:515–22.
Emerson BC. Alarm bells for the molecular clock? No support for Ho et al’.s model of time-dependent molecular rate estimates. Syst Biol. 2007;56:337–45.
Hsu F-H, Lin F-J, Lin Y-S. Phylogeographic structure of the Formosan wood mouse, Apodemus semotus Thomas. Zool Stud. 2001;40:91–102.
Michaux JR, Libois R, Paradis E, Filippucci M-G. Phylogeographic history of the yellow-necked fieldmouse (Apodemus flavicollis) in Europe and in the Near and Middle East. Mol Phylogenet Evol. 2004;32:788–98.
Michaux J, Bellinvia E, Lymberakis P. Taxonomy, evolutionary history and biogeography of the broad‐toothed field mouse (Apodemus mystacinus) in the eastern Mediterranean area based on mitochondrial and nuclear genes. Biol J Linn Soc. 2005;85:53–63.
Serizawa K, Suzuki H, Tsuchiya K. A phylogenetic view on species radiation in Apodemus inferred from variation of nuclear and mitochondrial genes. Biochem Genet. 2000;38:27–40.
Suzuki H, Sato JJ, Tsuchiya K, Luo J, Zhang Y, Wang Y, et al. Molecular phylogeny of wood mice (Apodemus, Muridae) in East Asia. Biol J Linn Soc. 2003;80:469–81.
Suzuki H, Yasuda SP, Sakaizumi M, Wakana S, Motokawa M, Tsuchiya K. Differential geographic patterns of mitochondrial DNA variation in two sympatric species of Japanese wood mice, Apodemus speciosus and A. argenteus. Genes Genet Syst. 2004;79:165–76.
Tomozawa M, Suzuki H. A trend of central versus peripheral structuring in mitochondrial and nuclear gene sequences of the Japanese wood mouse, Apodemus speciosus. Zool Sci. 2008;25:273–85.
Hayashi R, Takahara H, Tanida K, Danhara T. Vegetation response to East Asian monsoon fluctuations from the penultimate to last glacial period based on a terrestrial pollen record from the inland Kamiyoshi Basin, western Japan. Palaeogeogr Palaeoclimatol Palaeoecol. 2009;284:246–56.
Igarashi Y, Naruse T, Yatagai S, Danhara T. Vegetation and climate history since MIS 7 in the Kenbuchi Basin, northern Hokkaido :Comparison between MIS 6/5e and MIS 2/1. Quat Res (Japan) 2012;51:175–191. (in Japanese with English Summary).
Tomozawa M, Nunome M, Suzuki H, Ono H. Effect of founding events on coat colour polymorphism of Apodemus speciosus (Rodentia: Muridae) on the Izu Islands. Biol J Linn Soc. 2014;113:522–35.
Lisiecki LE, Raymo ME. A Pliocene‐Pleistocene stack of 57 globally distributed benthic δ18O records. Paleoceanography. 2005;20.
Sambrook J, Fritsch EF, Maniatis T. Molecular cloning: a laboratory manual. 2nd ed. New York: Cold Spring Harbor Laboratory Press; 1989.
Suzuki H, Nunome M, Kinoshita G, Aplin KP, Vogel P, Kryukov AP, et al. Evolutionary and dispersal history of Eurasian house mice Mus musculus clarified by more extensive geographic sampling of mitochondrial DNA. Heredity (Edinb). 2013;111:375–90.
Kocher TD, Thomas WK, Meyer A, Edwards SV, Pääbo S, Villablanca FX, et al. Dynamics of mitochondrial DNA evolution in animals: amplification and sequencing with conserved primers. Proc Natl Acad Sci. 1989;86:6196–200.
Yasuda SP, Vogel P, Tsuchiya K, Han S-H, Lin L-K, Suzuki H. Phylogeographic patterning of mtDNA in the widely distributed harvest mouse (Micromys minutus) suggests dramatic cycles of range contraction and expansion during the mid-to late Pleistocene. Can J Zool. 2005;83:1411–20.
Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011;28:2731–9.
Bandelt H-J, Forster P, Röhl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16:37–48.
Excoffier L, Lischer HEL. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10:564–7.
Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123:585–95.
Fu Y-X. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147:915–25.
Rogers AR, Harpending H. Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol. 1992;9:552–69.
Schneider S, Excoffier L. Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates vary among sites: application to human mitochondrial DNA. Genetics. 1999;152:1079–89.
Excoffier L. Patterns of DNA sequence diversity and genetic structure after a range expansion: lessons from the infinite‐island model. Mol Ecol. 2004;13:853–64.
Harpending HC. Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum Biol. 1994;66(4):591–600.
Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7:214.
Rambaut A, Drummond AJ. Tracer v1. 5. 2009. Available from http://beast.bio.ed.ac.uk.
Hewitt G. The genetic legacy of the Quaternary ice ages. Nature. 2000;405:907–13.
Nunome M, Torii H, Matsuki R, Kinoshita G, Suzuki H. The Influence of Pleistocene refugia on the evolutionary history of the Japanese hare, Lepus brachyurus. Zool Sci. 2010;27:746–54.
Ohnishi N, Uno R, Ishibashi Y, Tamate HB, Oi T. The influence of climatic oscillations during the Quaternary Era on the genetic structure of Asian black bears in Japan. Heredity (Edinb). 2009;102:579–89.
Nakagawa T, Kitagawa H, Yasuda Y, Tarasov PE, Nishida K, Gotanda K, et al. Asynchronous climate changes in the North Atlantic and Japan during the last termination. Science (80-). 2003;299:688–91.
Igarashi Y, Zharov AE. Climate and vegetation change during the late Pleistocene and early Holocene in Sakhalin and Hokkaido, northeast Asia. Quat Int. 2011;237:24–31.
Bintanja R, van de Wal RSW, Oerlemans J. Modelled atmospheric temperatures and global sea levels over the past million years. Nature. 2005;437:125–8.
Takahara H, Igarashi Y, Hayashi R, Kumon F, Liew P-M, Yamamoto M, et al. Millennial-scale variability in vegetation records from the East Asian Islands: Taiwan, Japan and Sakhalin. Quat Sci Rev. 2010;29:2900–17.
Tarasov P, Granoszewski W, Bezrukova E, Brewer S, Nita M, Abzaeva A, et al. Quantitative reconstruction of the last interglacial vegetation and climate based on the pollen record from Lake Baikal, Russia. Clim Dyn. 2005;25:625–37.
Tarasov PE, Nakagawa T, Demske D, Österle H, Igarashi Y, Kitagawa J, et al. Progress in the reconstruction of Quaternary climate dynamics in the Northwest Pacific: a new modern analogue reference dataset and its application to the 430-kyr pollen record from Lake Biwa. Earth-Science Rev. 2011;108:64–79.
Tarasov P, Bezrukova E, Karabanov E, Nakagawa T, Wagner M, Kulagina N, et al. Vegetation and climate dynamics during the Holocene and Eemian interglacials derived from Lake Baikal pollen records. Palaeogeogr Palaeoclimatol Palaeoecol. 2007;252:440–57.
Oshida T, Masuda R, Ikeda K. Phylogeography of the Japanese giant flying squirrel, Petaurista leucogenys (Rodentia: Sciuridae): implication of glacial refugia in an arboreal small mammal in the Japanese Islands. Biol J Linn Soc. 2009;98:47–60.
Oba T, Irino T. Sea level at the last glacial maximum, constrained by oxygen isotopic curves of planktonic foraminifera in the Japan Sea. J Quat Sci. 2012;27:941–7.
Goodman SJ, Tamate HB, Wilson R, Nagata J, Tatsuzawa S, Swanson GM, et al. Bottlenecks, drift and differentiation: the population structure and demographic history of sika deer (Cervus nippon) in the Japanese archipelago. Mol Ecol. 2001;10:1357–70.
Oshida T, Ikeda K, Yamada K, Masuda R. Phylogeography of the Japanese giant flying squirrel, Petaurista leucogenys, based on mitochondrial DNA control region sequences. Zool Sci. 2001;18:107–14.
Kawamoto Y, Shotake T, Nozawa K, Kawamoto S, Tomari K, Kawai S, et al. Postglacial population expansion of Japanese macaques (Macaca fuscata) inferred from mitochondrial DNA phylogeography. Primates. 2007;48:27–40.
Iwasa MA, Suzuki H. Evolutionary networks of maternal and paternal gene lineages in voles (Eothenomys) endemic to Japan. J Mammal. 2002;83:852–65.
Masson-Delmotte V, Hou S, Ekaykin A, Jouzel J, Aristarain A, Bernardo RT, et al. A review of Antarctic surface snow isotopic composition: observations, atmospheric circulation, and isotopic modeling*. J Clim. 2008;21:3359–87.
Ota Y. Coastal terraces of the Sado Island, Japan. Geogr Rev Jpn. 1964;37:226–42.
Yokoyama Y, Okuno J, Miyairi Y, Obrochta S, Demboya N, Makino Y, et al. Holocene sea‐level change and Antarctic melting history derived from geological observations and geophysical modeling along the Shimokita Peninsula, northern Japan. Geophys Res Lett. 2012;39:L13502.
Papadopoulou A, Anastasiou I, Vogler AP. Revisiting the insect mitochondrial molecular clock: the mid-Aegean trench calibration. Mol Biol Evol. 2010;27:1659–72.
Ho SYW, Phillips MJ, Cooper A, Drummond AJ. Time dependency of molecular rate estimates and systematic overestimation of recent divergence times. Mol Biol Evol. 2005;22:1561–8.
Burridge CP, Craw D, Fletcher D, Waters JM. Geological dates and molecular rates: fish DNA sheds light on time dependency. Mol Biol Evol. 2008;25:624–33.
Genner MJ, Seehausen O, Lunt DH, Joyce DA, Shaw PW, Carvalho GR, et al. Age of cichlids: new dates for ancient lake fish radiations. Mol Biol Evol. 2007;24:1269–82.
We are greatly indebted to Y. Igarashi, T. Irino, T. Iwasaki, G. Kinoshita, M. Motokawa, M. Sakaizumi, Y. Takada, and M. Yamamoto for providing valuable comments on an earlier version of the manuscript. We are also grateful for extensive comments from anonymous reviewers and the editor. This study was supported by a grant-in-aid for Scientific Research (B) to HS (no. 24405013) from the Japan Society for the Promotion of Science (JSPS).
The authors declare that they have no competing interests.
YS and HS designed the study; KT, YK, MT, and HS collected the data; and YS and MT performed the analyses. All authors contributed to writing the manuscript and approved its final version.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Suzuki, Y., Tomozawa, M., Koizumi, Y. et al. Estimating the molecular evolutionary rates of mitochondrial genes referring to Quaternary ice age events with inferred population expansions and dispersals in Japanese Apodemus . BMC Evol Biol 15, 187 (2015). https://doi.org/10.1186/s12862-015-0463-5
- Last Glacial Maximum
- Wood Mouse
- High Posterior Density
- Control Region Sequence
- Cytb Sequence