BACKGROUND
to investigate the recent evolutionary history of organisms, many phylogeographic studies have assessed molecular dating using rapidly evolving markers, such as mitochondrial dna  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 .
in fact, to analyze evolutionary history on a recent timescale, substantially higher rates  in the mtdna control region  have been proposed in studies of rodents such as the field vole  and house mouse , based on paleoclimatological and/or geological events using bayesian inferences, contrary to the rates “traditionally” used for mammals .
however, how long such a high evolutionary rate is applicable  <cit>  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  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 .
two endemic apodemus species in japan  are known to have a long  evolutionary history, independent from that of congeneric continental species according to molecular phylogenetic studies  <cit> .
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  <cit> .
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  was severe in the northern japanese islands, such as hokkaido and northern honshu  <cit> .
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  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  <cit> .
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  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  <cit> , 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.
methods
sample collection
in total, we used  <dig> a. argenteus individuals  and  <dig> a. speciosus individuals  .
the samples included dna and tissue samples stored in our laboratory  <cit> , in addition to new specimens from hokkaido, collected between  <dig> and  <dig>  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.fig.
<dig> collection locations for  apodemus argenteus and  a. speciosus.
the japanese islands were divided into five geographical regions, hokkaido, eastern honshu, western honshu, shikoku, and kyushu.
the solid line indicates the border of the karyotype forms of eastern  and western  honshu in a. speciosus.
circles indicate collection locations.
numbers with circles correspond to the locality number in table 1
apodemus argenteus
apodemus speciosus
<dig>  5321
b
sequences were obtained from the dna databases .
underlined dna codes are absent in the sequence of the control region
sequencing analysis
dna was extracted using a standard phenol-chloroform extraction method  <cit> .
for both species, two mtdna sequences were analyzed: the entire coding region of cytb  and half  of cr .
sequence determination for cytb was performed with semi-nested polymerase chain reaction , amplifying two fragments  in the second pcr  <cit> .
the first run was done using the primers l <dig> and h <dig>  <cit> , at a final concentration of  <dig>  pmol/μl, and amplitaq-gold dna polymerase  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  <dig>  pmol/μl for both upper  and lower  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 .
pcr for the cr sequences was performed with the amplitaq gold  <dig> master mix kit  using primers l <dig>  and cr <dig>  at a final concentration of  <dig>  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  and an abi <dig> automated sequencer .
the sequences of both strands were determined using the universal primers  in cytb and the primers used for pcr in cr.
the sequences were aligned using clustalw implemented in mega <dig>  <cit> .
phylogenetic analysis
phylogenetic trees with the concatenated mtdna sequences of cytb and cr were generated with mega <dig>  <cit> , using maximum-likelihood , maximum-parsimony , and neighbor-joining .
all indels were excluded from the phylogenetic analysis.
in addition to sequences of apodemus peninsulae , a. speciosus  and a. argenteus  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  was selected using mega <dig> with likelihood-ratio tests and the akaike information criterion .
the hasegawa–kishino–yano model  was used in the ml tree of a. speciosus.
in the nj analysis, the kimura two-parameter model  was used.
the reliability of nodes was assessed using  <dig> bootstrap replicates.
median-joining networks were constructed for both markers using the network program   <cit> .
we used arlequin  to calculate the nucleotide  and haplotype  diversity using all data  <cit> .
to detect rapid population expansion, mismatch distribution analyses, neutrality tests  were performed using arlequin .
the significance of neutrality was tested with  <dig>  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  <cit> .
the expected distribution was simulated under the sudden expansion model  <cit> .
we tested the validity of the sudden expansion model using a parametric bootstrap approach with  <dig>  replicates.
in this method, for each replicate, the sum of squared deviation  between the observed and expected distributions is compared with the ssd between the simulated and expected distributions using arlequin .
we used the raggedness index   <cit>  as a test statistic for the estimated sudden expansion models.
the expansion parameter tau  was estimated using arlequin   <cit>  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  and u is the cumulative evolutionary rate per generation for the whole sequence  <cit> .
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  for the hokkaido populations of both species was assumed to be 7– <dig> kyr bp, based on the timing of the recovery of the quercus population in northern hokkaido from the impact of the last glacial period  <cit> .
the time to the most recent common ancestor  of the cytb sequences and 95 % highest posterior density  were estimated using beast   <cit> .
the analysis was carried out with the hky <dig> substitution model and the strict-clock model using expansion growth as the tree prior.
two european species of apodemus sylvaticus  and apodemus flavicollis  whose divergence time was estimated at 2– <dig> million years ago  <cit> , were used as outgroup taxa.
bayesian searches were conducted using the markov chain monte carlo  method for  <dig> million generations.
the first  <dig> million generations were discarded as a burn-in.
independence among samplings was confirmed for each run by checking for high effective sample sizes .
tracer   <cit>  was used to assess convergence of the mcmc chains.
RESULTS
sequence data
sequence diversity indices for the two gene data sets are summarized in table  <dig>  complete cytb sequences  were obtained from  <dig> and  <dig> specimens of a. argenteus  and a. speciosus , respectively.
no deletion or insertion was detected in the cytb sequences of either species.
cr sequences for a. argenteus  and a. speciosus  contained  <dig> and  <dig> variable sites that defined  <dig> and  <dig> haplotypes, respectively.
in a. argenteus,  <dig> indels were detected, whereas only one insertion was detected in a. speciosus.table  <dig> genetic diversity indices of two mitochondrial dna markers  for two japanese apodemus species
s
h
hd
a. argenteus
a. speciosus
see figs.
<dig> and  <dig> for corresponding phylogroups of a. argenteus and a. speciosus, respectively
n sample size, s number of substitutions, h number of haplotypes, π nucleotide diversity, hd haplotype diversity
p <  <dig>  , p <  <dig>
population genetic analyses with mtdna
to determine the phylogeographic structures of the two apodemus species, phylogenetic trees and median-joining  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  revealed that the haplotypes fell into two clusters , which were supported by low  and high  bootstrap values depending the methods used, but they were unsupported  with the cr sequences alone .
the genetic distance for the cytb sequences between the two major clusters was calculated as  <dig>  % on average.
cluster i contained a relatively large subcluster, representing haplotypes from almost all japanese islands, from hokkaido to kyushu, with low bootstrap values .
the phylogenetic analyses further revealed a cluster representing all the hokkaido haplotypes  within subclade ia, with low to moderate supporting values .
notably, a haplotype from eastern honshu  clearly was closely associated with the hokkaido cluster ia- <dig>  although its supporting values were low .
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  and kyushu , and cluster iib included eight haplotypes from shikoku and kyushu.
the mj networks of the cytb, cr, and concatenated sequences  showed star-like structures for the haplotypes from hokkaido , as predicted in the phylogenetic tree analysis.
a portion of haplotypes from both western and eastern honshu also showed a star-like structure  in the cytb data set, although the structure was ambiguous in the networks of cr and concatenated data sets.fig.
<dig> maximum-likelihood trees for  apodemus argenteus and  a. speciosus from concatenated mitochondrial dna sequences .
numbers at nodes indicate bootstrap values by the maximum-likelihood, maximum-parsimonious, and neighbor-joining methods .
the taxon names correspond to those in table 1
fig.
<dig> median-joining networks constructed from  cytochrome b gene sequences  and  control region sequences , and  concatenated sequences  for  <dig> apodemus argenteus across japan.
the number of mutations  between haplotypes is indicated in branches.
the size of the circles is proportional to the number of samples.
the compositions of sample localities are also reflected in each haplotype.
the boxes with the solid and dashed lines are the major clusters  and clusters focused on in this study, respectively
the  <dig> haplotypes of the concatenated cytb and cr sequences from a. speciosus fell into two major clusters  that were always supported, but to varying degrees depending on the methods used .
cluster i included haplotypes from the three main islands of japan  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 , sado island , the izu islands , and the satsunan islands , 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  <dig>  %.
the mj networks of the cytb  and concatenated cytb and cr sequences  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 , 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 .
the mj network constructed from the concatenated sequences revealed two distinct groups in the haplotypes from hokkaido , with different levels of sequence divergence: iia- <dig>  and iia excluding iia- <dig> .
the iia- <dig> group was distributed throughout hokkaido, whereas the iiaex_iia- <dig> group was restricted to a few localities, such as rishiri island, kunashiri island, and the eastern part of hokkaido.
the sado island population  could be divided into subclusters iib- <dig> and iib- <dig>  the izu islands populations  into subclusters iic- <dig> and iic- <dig>  and the satsunan islands  populations into subclusters iid- <dig>  and iid- <dig> , each of which was further subdivided into two geographic areas.fig.
<dig> median-joining networks constructed from   <dig> cytochrome b gene sequences  and   <dig> control region sequences , and   <dig> concatenated sequences  for a. speciosus across japan.
the number of mutations  between haplotypes is indicated in branches.
the size of the circles is proportional to the number of samples.
the compositions of the sample localities are also reflected in each haplotype.
the boxes with solid and dashed lines are the major clusters  and the clusters addressed in this study, respectively
genetic diversity indices, including the nucleotide diversity  and haplotype diversity  values, were calculated for each of the phylogroups denoted in this study .
for a. argenteus, the π values of the total samples were  <dig>  % and  <dig>  %, and the hd values were  <dig>  and  <dig> , for cytb and cr, respectively.
notably, the π values of both markers for the hokkaido population were low .
for a. speciosus, the respective π values of the whole samples were  <dig>  % and  <dig>  %, and the hd values were  <dig>  and  <dig>  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 fs values for the hokkaido population of a. argenteus were significantly negative for both markers.
in a. speciosus, tajima’s d and fu’s fs values for the total sample set and the iia- <dig> group were significantly negative.
the iiaex_iia- <dig> had significant negative values for fu’s fs for both markers, although tajima’s d was not significant.
some other clusters, such as the ia- <dig> of a. argenteus  and cluster i of a. speciosus, also had significant negative values for the neutrality tests of both markers .
the mismatch distribution analysis of the cytb sequence was performed on each phylogroup except for two groups with small sample sizes : clusters i, ii, ia, ia- <dig>  and ia- <dig> for a. argenteus, and clusters i, ii, iia, iid, iia- <dig>  and iiaex_iia- <dig> for a. speciosus .
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- <dig> of cytb of a. speciosus .
in a. argenteus, analyses of the hokkaido population showed unimodal distributions with similar τ values for both markers .
the hokkaido population of a. speciosus tended to show bimodal distributions for both markers, although when separated, the two groups of hokkaido haplotypes, iia- <dig> and iiaex_iia- <dig>  had unimodal distributions for both markers, except cytb from iiaex_iia- <dig> fig.
<dig> mismatch distribution of the mitochondrial cytochrome b  and control region  sequences for a. argenteus  and a. speciosus .
bars indicate observed frequency, and a line denotes the expected frequency under the sudden expansion model.
ssd, sum of squared deviations; r, harpending’s raggedness index
ssd
r
τ
a. argenteus
cytb
a. speciosus
cytb
cytb
cytb
ssd sum of squared deviation, r harpending’s raggedness index, τ expansion parameter, t time  after expansion, μ evolutionary rate
*the t values for the predicted relatively recent expansion event of the hokkaido populations were set to  <dig>  and  <dig>  years ago based on the time frame of the recovery of the broadleaf forests in hokkaido  <cit> .
in addition, the possibility of expansion  immediately after the end of the lgm was also taken into account  <cit>
evolutionary rate estimation
the evolutionary rates  for both markers were calculated using the formula, t  = τ/2μk , suggesting that the hokkaido populations experienced rapid postglacial expansion, with the assumption that the latest expansion occurred 7– <dig> kyr bp, in association with the expansion of the quercus population in northern hokkaido  <cit> .
we also accounted for the possibility that the expansion of the wood mouse population occurred immediately after the lgm , given that apodemus includes habitat generalist species and that tree species in hokkaido largely recovered after the end of the glacial maximum  <cit> .
the evolutionary rates of cytb were estimated as  <dig> – <dig>  %/site/myr and  <dig> – <dig>  %/site/myr for the hokkaido populations of a. argenteus  and a. speciosus , respectively .
the evolutionary rates of the cr sequence were  <dig> – <dig>  %/site/myr and  <dig> – <dig>  %/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,  <dig> kyr bp, and obtained a value of  <dig>  %/site/myr .
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 , which suggested less regular substitution fixation rates over evolutionary time .
in the beast analysis for a. argenteus , the divergence time between the hokkaido population and the closest related haplotype from honshu  was estimated at  <dig> ,  <dig> , and  <dig>  kyr bp, using evolutionary rates of  <dig> ,  <dig>  and  <dig>  %/site/myr, respectively.
the other older divergences addressed in this study  were assessed using faster evolutionary rates of  <dig> ,  <dig> , and  <dig>  %/site/myr and the rather conservative evolutionary rates reported previously,  <dig>  %/site/myr  <cit>  and  <dig>  %/site/myr  <cit> , in addition to the value of  <dig>  %/site/myr.
the resulting estimated divergence times, including those of the two outgroup taxa, a. flavicollis and a. sylvaticus, are summarized in table  <dig> fig.
<dig>
a sea-level changes during the last  <dig> thousand years before present .
seven noticeable epochs of the sharp sea-level declines are marked with the stages of mis  <cit> .
b divergence time estimates of lineages representing the four peripheral insular domains of apodemus speciosus, based on a bayesian strict molecular clock applied to the mitochondrial cytochrome b  sequences .
the vertical bars indicate the 95 % highest posterior density intervals
a. argenteus
a. speciosus
a. flavicollis/a.
sylvaticus
four different evolutionary rates  were used.
tmrca and the 95 % highest posterior density  were computed with beast   <cit>  using the substitution model hky and the strict clock model.
the kimura 2-parameter distance  was calculated using mega <dig>  <cit> .
codes for clusters and nodes are as shown in figs.
<dig> and 5
discussion
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 .
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  and neutrality tests  suggested that both a. argenteus  and a. speciosus  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 .
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  <cit> .
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  <cit> .
the impact of the lgm and subsequent younger dryas glacial re-advancement  changed the tree species composition of hokkaido forests from cold coniferous species to temperate broad-leaved species, including oak trees   <cit> .
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 , although we can not exclude the idea that the expansion event occurred shortly after the lgm,  <dig> kyr bp.
the mismatch analysis and network analysis of hokkaido a. speciosus populations  suggested another expansion of iiaex_iia- <dig>  with a τ of  <dig> .
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 τ =  <dig>  is most likely to have occurred 7– <dig> kyr bp, the other expansion event for iiaex_iia- <dig> is likely attributable to the next nearest global climate fluctuation : more specifically, the periods  <dig> kyr bp or  <dig> kyr bp, equivalent to the terminations of marine isotope stages  <dig>  and  <dig> , respectively .
deciduous broad-leaved forests are known to have decreased and increased during these periods in japan  <cit> .
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  <dig>  as discussed later.
the haplotypes belonging to the iiaex_iia- <dig> 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 , 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  and mismatch distribution analysis  suggested a historical population expansion of a. speciosus in honshu, shikoku, kyushu, and some peripheral islands, such as the tsushima islands .
comparing the τ value with the hokkaido population , 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 .
pronounced changes in regional vegetation at the beginning of the last interglacial period  have been recorded in pollen-based quantitative biome reconstructions from the north eurasian study sites  <cit>  and central honshu  <cit> .
given that the postglacial appearance of quercus trees is evident  <cit> , the rapid expansion of a. speciosus  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  <cit> .
in addition, cluster i expansion generated the lineage of the tsushima islands ; this lineage is separated from kyushu by the tsushima strait, which is thought to have narrowed during the penultimate glacial maximum  <cit> .
thus, the expansion of cluster i conceivably started  <dig> 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 .
by contrast, except for the hokkaido subcluster, ia- <dig>  the network pattern of cluster i appears to have multiple subclusters , 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 —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  <cit> , black bears  <cit> , sika deer  <cit> , flying squirrels  <cit> , macaques  <cit> , and voles  <cit> .
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  <dig> –16 %/site/myr and  <dig> – <dig>  %/site/myr for a. argenteus and a. speciosus, respectively .
second, if we accept that the earlier expansion of the hokkaido population of a. speciosus started  <dig> kyr bp, the rate of cytb is  <dig>  %/site/myr.
finally, if we assume the expansion of a. specious in the southern islands of honshu, shikoku, and kyushu started at  <dig> kyr bp, the evolutionary rate is  <dig>  %/site/myr, which is comparable to those obtained previously  in the interspecies phylogenetic inferences with fossil evidence.
accordingly, these rate estimates can be categorized into faster  and slower  rates.
in a. speciosus, when we ran the beast analysis with the faster cytb evolutionary rate, the resulting estimates for divergences were quite recent .
for example, the divergences among the islands of izu, sado, and hokkaido  and between kyushu and the satsunan islands  were  <dig> and  <dig> kyr bp, respectively, suggesting dispersals in the mis  <dig> period .
however, this would seem highly unlikely because the sea-level drop at mis  <dig> is thought to have been at most around −90 m .
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  <dig> for interisland dispersals.
in contrast, the beast analysis with the lower evolutionary rate of  <dig>  %/site/myr  suggests that the historical dispersals of a. speciosus to the remote islands happened at or immediately after either mis  <dig> glacial periods  or older ones, of  <dig> kyr each  <cit> .
tmrcas for two lineages of izu , and sado , osumi , and tokara  islands were estimated to be  <dig> ,  <dig>  ,  <dig>  , and  <dig>  , respectively.
notably, the case of sado island is consistent with the geological evidence: at 110– <dig> kyr bp in the last interglacial period, sado island was separated into two islands by a sea strait  <cit> .
additionally, the resulting estimates of tmrcas of cluster i/ii, cluster ii, the lineages on hokkaido/sado/izu, and satsunan islands are  <dig> ,  <dig> ,  <dig> , and  <dig> , respectively .
these timings may have been associated with the glacial periods of mis  <dig> , mis  <dig> , and mis  <dig> , respectively.
the sea-level drops during these periods were more than 100 m, and thus, the peripheral 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  <dig> kyr bp ,  <dig> , and  <dig>  for clusters i and ii of a. argenteus from the slower evolutionary rates of  <dig>  %,  <dig>  %, and  <dig>  %, respectively .
thus, it seems reasonable to propose that the ancient divergence was associated with the end of glacial periods, most likely that following mis  <dig> .
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  was estimated as 22– <dig> kyr bp, with evolutionary rates of  <dig> –16 %/site/myr.
this indicates that the mtdna of a. argenteus came to hokkaido over the tsugaru strait  at most  <dig> kyr bp.
the sea level of the japan sea during the lgm is presumed to have been ~120 m lower than that at present  <cit> .
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 , 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  <cit> .
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 , whereas that calculated over a long timescale counts mutations that are fixed in the population by natural selection or genetic drift   <cit> .
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  <cit> .
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  <dig> kyr bp, is faster; e.g., 8–17 %/site/myr and 11–32 %/site/myr for cytb and cr, respectively .
in contrast, the use of the slower evolutionary rate  is more favorable for assessing older lineage divergences and population expansion events,  <dig> kyr bp or older.
this in turn indicates that the evolutionary rate declines dramatically within  <dig> kyr.
the presumed duration of the elevated evolutionary rate is shorter than those suggested for humans or insects   <cit> , and is comparable with that estimated in galaxiid fish   <cit>  and cichlid fish   <cit> .
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.
CONCLUSIONS
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 , which change depending on the time that has passed over the last  <dig>  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–ab <dig>  sequence data files in nexus file format, together with supplementary information files are available for downloading from the dryad repository: doi: <dig> /dryad.4902q.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
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.

 to investigate the recent evolutionary history of organisms, many phylogeographic studies have assessed molecular dating using rapidly evolving markers, such as mitochondrial dna 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. in fact, to analyze evolutionary history on a recent timescale, substantially higher rates in the mtdna control region have been proposed in studies of rodents such as the field vole and house mouse, based on paleoclimatological and/or geological events using bayesian inferences, contrary to the rates “traditionally” used for mammals. however, how long such a high evolutionary rate is applicable 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 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. two endemic apodemus species in japan are known to have a long evolutionary history, independent from that of congeneric continental species according to molecular phylogenetic studies. 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. 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 was severe in the northern japanese islands, such as hokkaido and northern honshu. 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 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. 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 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. methods sample collection in total, we used a. argenteus individuals and a. speciosus individuals. the samples included dna and tissue samples stored in our laboratory, in addition to new specimens from hokkaido, collected between and 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.fig. collection locations for apodemus argenteus and a. speciosus. the japanese islands were divided into five geographical regions, hokkaido, eastern honshu, western honshu, shikoku, and kyushu. the solid line indicates the border of the karyotype forms of eastern and western honshu in a. speciosus. circles indicate collection locations. numbers with circles correspond to the locality number in table 1 apodemus argenteus apodemus speciosus 5321 b sequences were obtained from the dna databases. underlined dna codes are absent in the sequence of the control region sequencing analysis dna was extracted using a standard phenol-chloroform extraction method. for both species, two mtdna sequences were analyzed: the entire coding region of cytb and half of cr. sequence determination for cytb was performed with semi-nested polymerase chain reaction, amplifying two fragments in the second pcr. the first run was done using the primers l and h, at a final concentration of pmol/μl, and amplitaq-gold dna polymerase 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 pmol/μl for both upper and lower 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. pcr for the cr sequences was performed with the amplitaq gold master mix kit using primers l and cr at a final concentration of 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 and an abi automated sequencer. the sequences of both strands were determined using the universal primers in cytb and the primers used for pcr in cr. the sequences were aligned using clustalw implemented in mega. phylogenetic analysis phylogenetic trees with the concatenated mtdna sequences of cytb and cr were generated with mega, using maximum-likelihood, maximum-parsimony, and neighbor-joining. all indels were excluded from the phylogenetic analysis. in addition to sequences of apodemus peninsulae, a. speciosus and a. argenteus 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 was selected using mega with likelihood-ratio tests and the akaike information criterion. the hasegawa–kishino–yano model was used in the ml tree of a. speciosus. in the nj analysis, the kimura two-parameter model was used. the reliability of nodes was assessed using bootstrap replicates. median-joining networks were constructed for both markers using the network program. we used arlequin to calculate the nucleotide and haplotype diversity using all data. to detect rapid population expansion, mismatch distribution analyses, neutrality tests were performed using arlequin. the significance of neutrality was tested with 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. we tested the validity of the sudden expansion model using a parametric bootstrap approach with replicates. in this method, for each replicate, the sum of squared deviation between the observed and expected distributions is compared with the ssd between the simulated and expected distributions using arlequin. we used the raggedness index as a test statistic for the estimated sudden expansion models. the expansion parameter tau was estimated using arlequin 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 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 for the hokkaido populations of both species was assumed to be 7– 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 of the cytb sequences and 95 % highest posterior density were estimated using beast. the analysis was carried out with the hky substitution model and the strict-clock model using expansion growth as the tree prior. two european species of apodemus sylvaticus and apodemus flavicollis whose divergence time was estimated at 2– million years ago, were used as outgroup taxa. bayesian searches were conducted using the markov chain monte carlo method for million generations. the first million generations were discarded as a burn-in. independence among samplings was confirmed for each run by checking for high effective sample sizes. tracer was used to assess convergence of the mcmc chains.  sequence data sequence diversity indices for the two gene data sets are summarized in table  complete cytb sequences were obtained from and specimens of a. argenteus and a. speciosus, respectively. no deletion or insertion was detected in the cytb sequences of either species. cr sequences for a. argenteus and a. speciosus contained and variable sites that defined and haplotypes, respectively. in a. argenteus, indels were detected, whereas only one insertion was detected in a. speciosus.table genetic diversity indices of two mitochondrial dna markers for two japanese apodemus species s h hd a. argenteus a. speciosus see figs.  and for corresponding phylogroups of a. argenteus and a. speciosus, respectively n sample size, s number of substitutions, h number of haplotypes, π nucleotide diversity, hd haplotype diversity p < , p <  population genetic analyses with mtdna to determine the phylogeographic structures of the two apodemus species, phylogenetic trees and median-joining 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 revealed that the haplotypes fell into two clusters, which were supported by low and high bootstrap values depending the methods used, but they were unsupported with the cr sequences alone. the genetic distance for the cytb sequences between the two major clusters was calculated as  % on average. cluster i contained a relatively large subcluster, representing haplotypes from almost all japanese islands, from hokkaido to kyushu, with low bootstrap values. the phylogenetic analyses further revealed a cluster representing all the hokkaido haplotypes within subclade ia, with low to moderate supporting values. notably, a haplotype from eastern honshu clearly was closely associated with the hokkaido cluster ia- although its supporting values were low. 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 and kyushu, and cluster iib included eight haplotypes from shikoku and kyushu. the mj networks of the cytb, cr, and concatenated sequences showed star-like structures for the haplotypes from hokkaido, as predicted in the phylogenetic tree analysis. a portion of haplotypes from both western and eastern honshu also showed a star-like structure in the cytb data set, although the structure was ambiguous in the networks of cr and concatenated data sets.fig. maximum-likelihood trees for apodemus argenteus and a. speciosus from concatenated mitochondrial dna sequences. numbers at nodes indicate bootstrap values by the maximum-likelihood, maximum-parsimonious, and neighbor-joining methods. the taxon names correspond to those in table 1 fig. median-joining networks constructed from cytochrome b gene sequences and control region sequences, and concatenated sequences for apodemus argenteus across japan. the number of mutations between haplotypes is indicated in branches. the size of the circles is proportional to the number of samples. the compositions of sample localities are also reflected in each haplotype. the boxes with the solid and dashed lines are the major clusters and clusters focused on in this study, respectively the haplotypes of the concatenated cytb and cr sequences from a. speciosus fell into two major clusters that were always supported, but to varying degrees depending on the methods used. cluster i included haplotypes from the three main islands of japan 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, sado island, the izu islands, and the satsunan islands, 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  %. the mj networks of the cytb and concatenated cytb and cr sequences 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, 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. the mj network constructed from the concatenated sequences revealed two distinct groups in the haplotypes from hokkaido, with different levels of sequence divergence: iia- and iia excluding iia-. the iia- group was distributed throughout hokkaido, whereas the iiaex_iia- group was restricted to a few localities, such as rishiri island, kunashiri island, and the eastern part of hokkaido. the sado island population could be divided into subclusters iib- and iib- the izu islands populations into subclusters iic- and iic- and the satsunan islands populations into subclusters iid- and iid-, each of which was further subdivided into two geographic areas.fig. median-joining networks constructed from cytochrome b gene sequences and control region sequences, and concatenated sequences for a. speciosus across japan. the number of mutations between haplotypes is indicated in branches. the size of the circles is proportional to the number of samples. the compositions of the sample localities are also reflected in each haplotype. the boxes with solid and dashed lines are the major clusters and the clusters addressed in this study, respectively genetic diversity indices, including the nucleotide diversity and haplotype diversity values, were calculated for each of the phylogroups denoted in this study. for a. argenteus, the π values of the total samples were  % and  %, and the hd values were and, for cytb and cr, respectively. notably, the π values of both markers for the hokkaido population were low. for a. speciosus, the respective π values of the whole samples were  % and  %, and the hd values were and 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 fs values for the hokkaido population of a. argenteus were significantly negative for both markers. in a. speciosus, tajima’s d and fu’s fs values for the total sample set and the iia- group were significantly negative. the iiaex_iia- had significant negative values for fu’s fs for both markers, although tajima’s d was not significant. some other clusters, such as the ia- of a. argenteus and cluster i of a. speciosus, also had significant negative values for the neutrality tests of both markers. the mismatch distribution analysis of the cytb sequence was performed on each phylogroup except for two groups with small sample sizes: clusters i, ii, ia, ia- and ia- for a. argenteus, and clusters i, ii, iia, iid, iia- and iiaex_iia- for a. speciosus. 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- of cytb of a. speciosus. in a. argenteus, analyses of the hokkaido population showed unimodal distributions with similar τ values for both markers. the hokkaido population of a. speciosus tended to show bimodal distributions for both markers, although when separated, the two groups of hokkaido haplotypes, iia- and iiaex_iia- had unimodal distributions for both markers, except cytb from iiaex_iia- fig. mismatch distribution of the mitochondrial cytochrome b and control region sequences for a. argenteus and a. speciosus. bars indicate observed frequency, and a line denotes the expected frequency under the sudden expansion model. ssd, sum of squared deviations; r, harpending’s raggedness index ssd r τ a. argenteus cytb a. speciosus cytb cytb cytb ssd sum of squared deviation, r harpending’s raggedness index, τ expansion parameter, t time after expansion, μ evolutionary rate *the t values for the predicted relatively recent expansion event of the hokkaido populations were set to and  years ago based on the time frame of the recovery of the broadleaf forests in hokkaido. in addition, the possibility of expansion immediately after the end of the lgm was also taken into account evolutionary rate estimation the evolutionary rates for both markers were calculated using the formula, t  = τ/2μk, suggesting that the hokkaido populations experienced rapid postglacial expansion, with the assumption that the latest expansion occurred 7– 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, 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 –  %/site/myr and –  %/site/myr for the hokkaido populations of a. argenteus and a. speciosus, respectively. the evolutionary rates of the cr sequence were –  %/site/myr and –  %/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, kyr bp, and obtained a value of  %/site/myr. 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, which suggested less regular substitution fixation rates over evolutionary time. in the beast analysis for a. argenteus, the divergence time between the hokkaido population and the closest related haplotype from honshu was estimated at,, and kyr bp, using evolutionary rates of, and  %/site/myr, respectively. the other older divergences addressed in this study were assessed using faster evolutionary rates of,, and  %/site/myr and the rather conservative evolutionary rates reported previously,  %/site/myr and  %/site/myr, in addition to the value of  %/site/myr. the resulting estimated divergence times, including those of the two outgroup taxa, a. flavicollis and a. sylvaticus, are summarized in table  fig. a sea-level changes during the last thousand years before present. seven noticeable epochs of the sharp sea-level declines are marked with the stages of mis. b divergence time estimates of lineages representing the four peripheral insular domains of apodemus speciosus, based on a bayesian strict molecular clock applied to the mitochondrial cytochrome b sequences. the vertical bars indicate the 95 % highest posterior density intervals a. argenteus a. speciosus a. flavicollis/a. sylvaticus four different evolutionary rates were used. tmrca and the 95 % highest posterior density were computed with beast using the substitution model hky and the strict clock model. the kimura 2-parameter distance was calculated using mega. codes for clusters and nodes are as shown in figs.  and 5 discussion 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. 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 and neutrality tests suggested that both a. argenteus and a. speciosus 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. 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. the impact of the lgm and subsequent younger dryas glacial re-advancement changed the tree species composition of hokkaido forests from cold coniferous species to temperate broad-leaved species, including oak trees. 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, although we can not exclude the idea that the expansion event occurred shortly after the lgm, kyr bp. the mismatch analysis and network analysis of hokkaido a. speciosus populations suggested another expansion of iiaex_iia- with a τ of. 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 τ =  is most likely to have occurred 7– kyr bp, the other expansion event for iiaex_iia- is likely attributable to the next nearest global climate fluctuation: more specifically, the periods kyr bp or kyr bp, equivalent to the terminations of marine isotope stages and, respectively. deciduous broad-leaved forests are known to have decreased and increased during these periods in japan. 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 as discussed later. the haplotypes belonging to the iiaex_iia- 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, 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 and mismatch distribution analysis suggested a historical population expansion of a. speciosus in honshu, shikoku, kyushu, and some peripheral islands, such as the tsushima islands. comparing the τ value with the hokkaido population, 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. pronounced changes in regional vegetation at the beginning of the last interglacial period have been recorded in pollen-based quantitative biome reconstructions from the north eurasian study sites and central honshu. given that the postglacial appearance of quercus trees is evident, the rapid expansion of a. speciosus 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; 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 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. by contrast, except for the hokkaido subcluster, ia- the network pattern of cluster i appears to have multiple subclusters, 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 —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 –16 %/site/myr and –  %/site/myr for a. argenteus and a. speciosus, respectively. second, if we accept that the earlier expansion of the hokkaido population of a. speciosus started kyr bp, the rate of cytb is  %/site/myr. finally, if we assume the expansion of a. specious in the southern islands of honshu, shikoku, and kyushu started at kyr bp, the evolutionary rate is  %/site/myr, which is comparable to those obtained previously in the interspecies phylogenetic inferences with fossil evidence. accordingly, these rate estimates can be categorized into faster and slower rates. in a. speciosus, when we ran the beast analysis with the faster cytb evolutionary rate, the resulting estimates for divergences were quite recent. for example, the divergences among the islands of izu, sado, and hokkaido and between kyushu and the satsunan islands were and kyr bp, respectively, suggesting dispersals in the mis period. however, this would seem highly unlikely because the sea-level drop at mis is thought to have been at most around −90 m. 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 for interisland dispersals. in contrast, the beast analysis with the lower evolutionary rate of  %/site/myr suggests that the historical dispersals of a. speciosus to the remote islands happened at or immediately after either mis glacial periods or older ones, of kyr each. tmrcas for two lineages of izu, and sado, osumi, and tokara islands were estimated to be,,, and, respectively. notably, the case of sado island is consistent with the geological evidence: at 110– 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,,, and, respectively. these timings may have been associated with the glacial periods of mis, mis, and mis, respectively. the sea-level drops during these periods were more than 100 m, and thus, the peripheral 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 kyr bp,, and for clusters i and ii of a. argenteus from the slower evolutionary rates of  %,  %, and  %, respectively. thus, it seems reasonable to propose that the ancient divergence was associated with the end of glacial periods, most likely that following mis. 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 was estimated as 22– kyr bp, with evolutionary rates of –16 %/site/myr. this indicates that the mtdna of a. argenteus came to hokkaido over the tsugaru strait at most 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, 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. 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, whereas that calculated over a long timescale counts mutations that are fixed in the population by natural selection or genetic drift. 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 kyr bp, is faster; e.g., 8–17 %/site/myr and 11–32 %/site/myr for cytb and cr, respectively. in contrast, the use of the slower evolutionary rate is more favorable for assessing older lineage divergences and population expansion events, kyr bp or older. this in turn indicates that the evolutionary rate declines dramatically within kyr. the presumed duration of the elevated evolutionary rate is shorter than those suggested for humans or insects, and is comparable with that estimated in galaxiid fish and cichlid fish. 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, which change depending on the time that has passed over the last  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–ab sequence data files in nexus file format, together with supplementary information files are available for downloading from the dryad repository: doi: /dryad.4902q. competing interests the authors declare that they have no competing interests. authors’ contributions 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. 
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. assuming that global warming after the last glacial period 7– thousand years before present was associated with the expansion, the evolutionary rates 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. 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. we examined mitochondrial cytochrome b and control region 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. 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  %/site/myr under the assumption of a postglacial population expansion event long ago, most probably at kyr bp. however, the slower evolutionary rate of  %/site/myr would be reasonable for several divergence events that were associated with glacial periods older than 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. the elevated evolutionary rate seemed to decline within  years.
 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 and control region 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– thousand years before present was associated with the expansion, the evolutionary rates 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  %/site/myr under the assumption of a postglacial population expansion event long ago, most probably at 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  %/site/myr would be reasonable for several divergence events that were associated with glacial periods older than 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  years. issue-copyright-statement© the author 2015 