Late Quaternary Environmental and Human Impacts on the Mitochondrial DNA Diversity of Four Commensal Rodents in Myanmar

We addressed the spatiotemporal characteristics of four commensal rodent species occurring in Myanmar in comparison with other areas of the Indo-Malayan region. We examined sequence variations of the mitochondrial cytochrome b gene ( Cytb ) in the Pacific rat ( Rattus exulans ), roof rat ( Rattus rattus complex, RrC), lesser bandicoot rat ( Bandicota bengalensis ), and house mouse ( Mus musculus ) using the recently developed time-dependent evolutionary rates of mtDNA. The Cytb sequences of RrC from Myanmar were shown to belong to RrC Lineage II, and their level of genetic diversity was relatively high compared to those of the other three species. RrC was found to have experienced bottleneck and rapid expansion events at least twice in the late Pleistocene period in Myanmar and a nearby region. Accordingly, paleoclimatic environmental fluctuations were shown to be an important factor affecting rodents in the subtropics of the Indo-Malayan region. Our results show that human activities during the last 10,000 years of the Holocene period affected the population dynamics of the rodent species examined, including introducing them to Myanmar from neighboring countries. Further study of these four commensal rodents in other geographic areas of the Indo-Malayan region would allow us to better understand the factors that drove their evolution and their ecological trends.


Introduction
The subfamily Murinae is a large taxonomic group that consists of more than 500 species found mainly in the Indo-Malayan region, a region of southern Asia with a high level of terrestrial animal biodiversity (Musser and Carleton 2005;Srinivasulu and Srinivasulu 2012) covering South and Southeast Asia and southern China. Some murine rodents, such as rats and mice, that live within or close to human habitations are called commensal rodents and have affected human activities in a variety of ways, such as by causing agricultural damage (Aplin et al. 2003) and carrying infectious diseases (Kosoy et al. 2015). On the other hand, such species can provide useful information for clarifying human movements from prehistoric to modern times (Boursot et al. 1993;Matisoo-Smith and Robins 2004;Aplin et al. 2011;Jones et al. 2013; Thomson et al. 2014). Phylogeographic analyses have been conducted for a number of species of commensal rodents, and anthropological factors were found to have shaped the genetic structures of commensal rodent species over the last 15,000 years (Auffray and Britton-Davidian Electronic supplementary material The online version of this article (https://doi.org/10.1007/s10914-020-09519-4) contains supplementary material, which is available to authorized users. 2012; Aplin et al. 2011;Conroy et al. 2013;Suzuki et al. 2013). The impacts of both humans and the natural environment were likely involved in the processes that shaped their genetic variation. However, the evolutionary histories of many of these species have not yet been fully elucidated.
Myanmar is located in the central part of the Indo-Malayan region and plays a role connecting neighboring countries including India, China, Thailand, and Indonesia. In addition, Myanmar harbors its own species and lineages, and it is essential to understand the functions of physical barriers such as mountain ranges in the emergence of genetic variation within these lineages. In the Indo-Malayan region, the Arakan Mountain system on the border of India and Myanmar may act as a geographic barrier, dividing the region into two parts (Chingangbam et al. 2015). Mountain ranges have been shown to act as borders between species or local intraspecies phylogroups in species of the genus Mus (Shimada et al. 2007(Shimada et al. , 2009(Shimada et al. , 2010Myat Myat Zaw et al. 2019). Widespread taxonomic groups are ideal systems for testing biogeographical hypotheses, as they provide an opportunity for examining how geographical barriers, as well as the processes of dispersal and colonization, affect diversification and determine biogeographical distribution patterns (e.g., Irschick et al. 1997). In Myanmar, commensal murine rodents are abundant and represented by the Pacific rat (Rattus exulans), roof rat (Rattus rattus complex, RrC), lesser bandicoot rat (Bandicota bengalensis), and house mouse (Mus musculus). Due to their large impacts on human welfare (Thomson et al. 2018), it is necessary to elucidate the genetic structure and inferred demographics of commensal rodents in Myanmar.
Rattus exulans is a common commensal rodent variably described as a 'ground rat' (Medway 1969) and 'highly arboreal' (Aplin et al. 2003). It is distributed in Southeast Asia, the Indonesian islands, and the Pacific Islands. According to Corbet and Hill (1992), it is a common pest species found in houses, granaries, cultivated areas, scrubland, and forests, usually on the ground. In Thailand, it has been reported in houses, markets, the vicinity of houses, and in rice paddies located near buildings (Marshall 1988). Rattus exulans has also attracted the attention of researchers interested in anthropogenic movements across the deep sea via ships in the Pacific Islands (Matisoo-Smith and Robins 2004), with this species acting as a symbiotic animal by providing nutrition to humans. Human-mediated dispersal of these rats among the Pacific Islands is an issue of great concern. Based on a mitochondrial DNA (mtDNA) phylogeographic survey of R. exulans diversity across its potential natural range of mainland (MSEA) and island (ISEA) areas of Southeast Asia, prehistoric MSEA and ISEA rodent fauna have been confirmed (Thomson et al. 2014). The highest mtDNA lineage diversities, including significant haplotype diversity, were observed in ISEA populations of R. exulans. From the results of that study and other evidence, Flores Island, Indonesia, is predicted to be the origin of this species (Thomson et al. 2014). The survey of MSEA remains limited to date, and it is necessary to address this issue through the collection of more specimens from MSEA.
The roof rat (or ship rat) exhibits strong arboreal features and is predicted to have its origins in Asia; thus, it is common in the Indo-Malayan region today. The species is known to include complicated geographic groups and was recently recognized as the RrC, containing six distinct mitochondrial cytochrome b (Cytb) sequences (Lineage I-VI) with apparent geographic affinity (Chinen et al. 2005;Pagès et al. 2010;Tollenaere et al. 2010;Aplin et al. 2011). Eastern and Southern India, representing RrC Lineage I, is the likely geographical focus of the 2n = 38 karyotype group, whereas the geographic area covering the easternmost part of India to southern China represents the natural range of RrC Lineage II, representing chromosome form 2n = 42, also called R. tanezumi. Although this mtDNA lineage extends to the easternmost part of Eurasia-the Japanese Islands (Chinen et al. 2005;Kambe et al. 2012), the place of origin of this species and evolutionary events of range expansion of RrC II have not yet been determined. Although phylogenetic and phylogeographic studies have been carried out for RrC (Pagès et al. 2010(Pagès et al. , 2013Aplin et al. 2011), intensive genetic analysis of populations in Myanmar has not yet been reported.
The bandicoot rats (genus Bandicota) are pest rodent species that are abundant and widespread throughout the continental portion of the Indo-Malayan region. Recently, a molecular phylogenetic study was conducted on the three bandicoot rat species from Myanmar, where all three species of this genus, B. bengalensis, B. indica, and B. savilei, occur and a number of genetic sequences of B. bengalensis from Myanmar have been published (Mori et al. 2020). While B. bengalensis inhabits a wide area of Myanmar and is relatively dense in population, the recovered Cytb sequences were shown to contain limited nucleotide diversity and to exhibit a localityspecific pattern; i.e., each locality tends to possess a single haplotype (Mori et al. 2020). However, no detailed temporal mechanism has been inferred to assess the evolutionary history of this species.
The evolutionary history of M. musculus has been extensively studied. Three major subspecies groups, including M. musculus castaneus (CAS), are known to have extended their ranges alongside prehistoric human movements during the development of agriculture across the Eurasian Continent (e.g., Boursot et al. 1993). The South Asian group M. m. castaneus is known to contain four distinct mtDNA sublineages, referred to here as CAS 1-4 (Suzuki et al. 2013). One of these lineages, CAS-1, is now present over a wide area of India, Southeast Asia, East Asia, and the associated islands of Japan, Philippines, and Indonesia after an episodic expansion event ca. 7000-8000 years ago (Suzuki et al. 2013;Kuwayama et al. 2017). The CAS-1 sublineage experienced a secondary expansion event in southern China that extended its range to Japan and southern Sakhalin (Suzuki et al. 2013). In a previous study, the extent of genetic diversity in Cytb sequences was found to be relatively small in Myanmar (Suzuki et al. 2013;Myat Myat Zaw et al. 2019), indicating that recent structuring of the population had occurred due to the introduction of human agricultural systems in Southeast Asia including Myanmar. However, details of the evolutionary progress of M. musculus in Myanmar are not yet known.
The influence of Quaternary environmental changes accompanying glacial-interglacial cycles and sea-level fluctuations (Lambeck et al. 2002;Lisiecki and Raymo 2005;Bintanja et al. 2005;Jansen et al., 2017), including major impacts shaping the genetic diversity of terrestrial animals, have been well documented. In particular, the last (LGM, ca. 20,000 years ago) and penultimate (PGM, ca. 140,000 years ago) glacial maxima are reported to have had strong effects on the present genetic structures of species with temperate origins. Population expansion events have been deduced for rodents dwelling in temperate forests (genus Apodemus) (Hanazaki et al. 2017). In addition, population expansion events have been inferred for subarctic grassland-dwelling rodents (genus Myodes) (Kohli et al. 2015;Honda et al. 2019) during the intermediate period 50,000-60,000 years ago, when warmer conditions were beginning during early marine oxygen isotope stage 3 (MIS 3) after a substantially colder period (MIS 4) (e.g., Sawagaki et al. 2004;Shang et al. 2018). Therefore, we sought to determine whether these critical periods of drastic environmental change during the late Quaternary affected the population dynamics of species in subtropical areas.
In this study, we addressed the evolutionary history of four commensal rodents that are common in Myanmar based on variations in their mtDNA sequences. We collected R. exulans, RrC, and M. musculus from Myanmar and analyzed their Cytb nucleotide sequences. We also obtained Cytb sequences for R. exulans (Thomson et al. 2014), M. musculus (e.g., Myat Myat Zaw et al. 2019), and B. bengalensis (Mori et al. 2020) and sequences for RrC and M. musculus from Manipur, a state in northeastern India (Chingangbam et al. 2015) from databases. In addition, Cytb sequences from other geographic regions were obtained from databases for phylogeographic assessment. The study of these four commensal rodent species allows us to compare their histories and discuss their commonalities and differences. This in turn may provide useful clues to better understand the factors driving genetic variation and to elucidate the zoogeographic role of Myanmar, the center of the Indo-Malayan region.

Sample Collection
Sample collection in Myanmar was conducted mainly by staff members of the Department of Zoology, University of Yangon from 2013 to 2018. For this study, 35 individuals of Rattus exulans, 13 individuals of RrC, and 4 individuals of M. musculus were collected and subjected to molecular analysis (Table 1, Fig. 1).

Sequencing Analysis
DNA was extracted from the liver tissue from each individual using a Qiagen DNA Isolation kit (Hilden, Germany) following the manufacturer's protocol. Sequence determination was performed for Cytb (1034-1140 bp) in accordance with a method described previously . The gene region was amplified through polymerase chain reaction (PCR) using the following protocol: 95°C for 10 min, followed by 30 cycles of 95°C for 30 s, 50°C for 30 s, and 60°C for 60 s, and finally, 72°C for 7 min. PCR products were seq u e n c e d u s i n g a P R I S M R e a d y R e a c t i o n D y e DeoxyTerminator Cycle Sequencing kit Ver 3.1 and an ABI 3500 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA). The nucleotide sequences reported in the present study appear in the DDBJ, EMBL, and GenBank nucleotide sequence databases under accession numbers LC510767-LC510818.

Phylogenetic Analysis
Long (1034 bp) and short (381 bp) Cytb sequences of R. exulans were obtained from databases and used for phylogenetic inference based on datasets of the MSEA and MSEA plus ISEA haplotypes. Database sequences of RrC (1140 bp, 1034 bp) were used in this study. We obtained 24 Cytb sequences of B. bengalensis (1026 bp, 1140 bp) from the databases, 23 of which were collected from Myanmar (MK654761-MK654784) and one from Pakistan (AM408336). Cytb sequences of M. musculus from Myanmar were obtained from databases (LC426295-LC426303). Phylogenetic trees were constructed based on the maximum likelihood (ML) method using MEGA ver. 7.0 software (Kumar et al. 2 0 1 6 ) w i t h t h e s u b s t i t u t i o n m o d e l s G T R + G ( R . e x u l a n s ) , T N 9 3 + G + I ( R r C ) , T N 9 3 + G (B. bengalensis), and TN93 (M. musculus). The best fit model was determined using the Akaike information criterion (AIC; Posada and Buckley 2004), as implemented in MEGA. Levels of genetic variation and divergence were assessed based on nucleotide diversity (π) using ARLEQUIN ver. 3.5 software (Excoffier and Lischer 2010).

Population genetic analysis
Colonization (founder effect) and post-glacial environmental impacts were assessed by detecting the rapid expansion signal (Excoffier and Schneider 1999). Networks were constructed using the Median-Joining (MJ) method (Bandelt et al. 1999), as implemented in PopART software (Leigh and Bryant 2015). Mismatch distribution analysis and the neutrality tests Tajima's D and Fu's F S were performed using ARLEQUIN (Excoffier and Lischer 2010). The expected distribution was simulated for the sudden expansion model through mismatch distribution analysis, in which the validity of the model was tested using a parametric bootstrap approach with 1000 replicates. In this method, the sum of squared deviation (SSD) between the observed distribution and the expected distribution was compared with the SSD between the simulated distribution and expected distribution for each replicate. The raggedness index (r; Harpending 1994) was used as a test statistic for the predicted sudden expansion model.
The expansion parameter tau (τ) was estimated for the Cytb datasets using ARLEQUIN (Excoffier and Lischer 2010) in clusters that showed signs of sudden demographic expansion. The 95% confidence interval (CI) around the estimated τ parameter was determined using the bootstrap approach. The time since expansion in generation (t) can be estimated by the formula t = τ/ 2uk, where u is the evolutionary rate per generation and k is sequence length (Rogers and Harpending 1992;Rogers 1995). Here, the time since expansion in years (T) was estimated by the formula T = τ/2μk, where μ is the evolutionary rate per year, since u = μg, and t = T/g, where g is generation time in years. Here, we used evolutionary rates (μ') per million years (myr), for convenience.
In this study, we examined the evolutionary rate of the mitochondrial gene sequence change in a timedependent manner (Ho et al. 2011). Previously, we addressed rapid expansion events seen in local populations of wood mice and voles occurring in the Japanese archipelago and applied biogeographic events as

Rattus exulans
For R. exulans, we compared the Cytb sequences (1034 bp) of 35 rats from Myanmar and 11 rats from other countries. A close relationship between haplotypes from Myanmar and those from other countries was observed in the ML tree (Fig. 2a) and in the MJ network (Fig. 3a), consistent with the low level of nucleotide diversity (π = 0.0018; Table 2). Although the pattern obtained for the MJ network exhibited a broken star-like structure, the  Fig. S1) and neutrality tests did not support the rapid expansion model (Table 2). Phylogeographic analysis employed a 381-bp Cytb sequence dataset using database sequences (e.g., Thomson et al. 2014) covering haplotypes from ISEA (n = 10) and MSEA (n = 24). A star-like network was obtained, similar to that reported previously in a study of the mtDNA control region (Fig. 3a; Thomson et al. 2014). The neutrality tests provided evidence of rapid expansion. The mismatch distribution analysis resulted in a τ value of 0.94 (CI: 0.07-1.97) ( Table 2). The time of expansion initiation was estimated as 11,200 years ago (CI: 8400-23,500 years ago) using an evolutionary rate of 0.11 substitutions/site/myr.

RrC
We collected 13 rats from Kalaymyo, Nay Pyi Taw, Pyay, Yangon, and Pathein along the Ayeyarwady River and found that all rats from Myanmar had pure white hairs on the ventral side. Analysis of the Cytb sequences (1140 bp) yielded eight haplotypes in total. ML tree construction of these and available database sequences (1034 bp) revealed that the haplotypes from Myanmar all belonged to RrC Lineage II (Fig. 2b), one of the six major mtDNA lineages of RrC defined by Aplin et al. (2011). This lineage is sometimes called 'R. tanezumi' and represents rats from East Asia, including the Japanese Islands (Fig. 1b).
We constructed an MJ network using Cytb sequences (1140 bp) from Myanmar and other countries (n = 30; Fig.  3b). Network construction revealed a cluster with a star-like

Bandicota bengalensis
We obtained 24 Cytb sequences (1026 bp) for B. bengalensis from the databases, consisting of 23 from Myanmar (Mori et al. 2020) and one from Pakistan. We then constructed an ML tree (Fig. 2c) and MJ network, in which an apparent starlike structure was not seen (Fig. 3c). Mismatch distribution analysis ( Supplementary Fig. S1) and neutrality tests were performed on the 23 sequences from Myanmar; the Tajima's D obtained was significantly negative (P = 0.014), and Fu's Fs was not significant ( Table 2). The τ value was calculated to be 0.91 (CI: 0.43-1.52), with an estimated time of expansion initiation, using the evolutionary rate of 0.11 substitutions/ site/myr, of 4000 (CI: 1900-6700) years ago.

Mus musculus
We obtained four new Cytb sequences (1140 bp) of M. musculus from Taunggyi, Myanmar and combined those with previously determined sequences including three from Manipur, India (n = 19). We constructed an ML tree using these (Fig. 2d) and MJ network (Fig. 3d), revealing that all 15 haplotypes from Myanmar and Manipur belonged to cluster CAS-1, the major lineage found over a large area of Southeast Asia and East Asia, including the Japanese Islands (Suzuki et al. 2013, Kuwayama et al. 2017. We assessed the population dynamics of M. musculus and obtained consistent signals of population expansion in the analyses of mismatch distribution ( Supplementary  Fig. S1) and neutrality tests ( Table 2). The τ values calculated for two datasets, Myanmar and Manipur (n = 19) and Myanmar plus other localities (n = 23), were 1.25 (CI: 0.12-1.97) and 1.36 (CI: 0.57-2.01), respectively. The initiation dates of the expansion events were calculated as 5000 (CI: 500-7900) and 5400 (CI: 2300-8000) years ago, respectively (Table 2).

Discussion
We conducted this study to elucidate the evolutionary histories of the Myanmar populations of four species of rodents that live in close proximity to humans-R. exulans, RrC, B. bengalensis, and M. musculus. By comparing patterns in their genetic structure and population dynamics, we noticed apparent differences among these four species in terms of their levels of genetic diversity and presence or absence of rapid expansion signals. Our data allowed us to reconstruct demographic histories and estimate times at which the expansion events began, although the accuracy of these calculations was limited by relatively small sample sizes. RrC exhibited relatively high genetic diversity, whereas the remaining three species possessed relatively low levels of genetic diversity (Fig.  3). Moreover, signals of rapid expansion events were observed for three species, RrC, B. bengalensis, and M. musculus, whereas no expansion signal was observed for R. exulans, using the dataset of haplotypes from Myanmar.
Our results allowed us to reconstruct their natural histories and to better understand the factors shaping the diversity of commensal rodents, which are assumed to operate in a speciesspecific manner (Fig. 4), as discussed below. The Myanmar population of R. exulans is characterized by a low level of genetic diversity, no expansion signal, no apparent geographic differentiation in Myanmar, and substantial affinity with the populations of neighboring countries. These results indicate that the population of R. exulans emerged relatively recently in Myanmar. Thomson et al. (2014) first revealed the deeper divergence of R. exulans populations in ISEA compared to those in MSEA, hypothesizing that this species originated in the Indonesian islands. The network that we created based on the dataset of short sequences (381 bp, Thomson et al. 2014) exhibited a phylogeographic trend (Fig.  3a) that supports the hypothesis of Thomson et al. (2014), exhibiting a rapid expansion signal, as has been reported previously (Thomson et al. 2014), and an expansion start time of ca. 11,000 years ago (Table 2, Fig. 4). The point in time when R. exulans arrived in Myanmar from other geographic regions such as Indonesia remains unclear. One possibility is that some haplotypes of ISEA origin were introduced quite recently, even in modern times.
For RrC, the mtDNA sequences we obtained from Myanmar were found to belong to RrC Lineage II, which is assumed to be distributed from Myanmar in the west to the Japanese Islands in the east (Aplin et al. 2011). We disclosed a distinct sublineage, referred to here as IIa-1, comprising haplotypes from Myanmar (this study) and the neighboring region of Manipur, India (Chingangbam et al. 2015). Hence, it has been revealed that the Arakan Mountain system on the border of India and Myanmar may not serve as a geographic barrier for roof rats. The cluster exhibited a star-like structure and possessed significant signals of rapid expansion, with a τ value of 6.79 (1140 bp; Table 2). Using an evolutionary rate of 0.047 substitutions/site/myr for Cytb sequences (Honda et al. 2019), the expansion time was estimated as 63,000 years BP (CI: 44,000-77,000 years ago), corresponding to the rapid warming period during early MIS 3, following a bottleneck event during the colder period of MIS 4 (e.g., Woltering et al. 2011) (Fig. 4). Notably, the entire set of haplotypes from Myanmar and Manipur exhibited expansion signals (Table 2), with a τ value of 8.06 and an estimated expansion start time of ca. 140,000 years ago, corresponding to MIS 5e. The cluster of haplotypes IIa (Fig. 2a, "IIA" in Aplin et al. 2011) from the easternmost countries of Eurasia, Vietnam, China, and Japan also exhibited evidence of expansion, with a τ value of 1.96 (Aplin et al. 2011; Table 2). These results further indicate that RrC Lineage II experienced a bottleneck and subsequent population expansion at three different time points, including during the influential periods of MIS 5e and early MIS 3 in the region of the Arakan Mountains (Fig. 4). The sudden expansion of IIa-2 in the eastern parts of Eurasia occurred as recently as 9400 years ago, the beginning of Holocene (Table 2), but it remains uncertain whether this was a result of human activities (Aplin et al. 2011) or of global environmental fluctuation.
Similar to R. exulans, B. bengalensis exhibits a low level of genetic diversity, as previously reported (Mori et al. 2020;Fig. 3, Table 2), although the species is quite common in paddy fields and livestock farms in Myanmar. Analogous to the above mentioned predicted history of R. exulans, the Fig. 4 Schematic representation of possible population expansion and dispersal processes inferred from mtDNA sequence variation in the Myanmar populations of Rattus exulans, Rattus rattus complex (RrC), Bandicota bengalensis, and Mus musculus. a The marine oxygen isotope curve for the last 140,000 years is adapted from Petit et al. (1999). b Rapid expansion events inferred from mtDNA analysis are plotted at the corresponding Marine Isotope Stage (MIS) periods: MIS 5e (RrC sequences from Myanmar and Manipur, India), early MIS 3 (RrC Lineage IIa-1), and early MIS 1 (RrC sequences from Vietnam, China, and Japan; see Aplin et al. 2011). The expansion events suggested by the Cytb sequences of R. exulans (Myanmar), B. bengalensis (Myanmar), and M. musculus (Myanmar and Manipur) are indicated. The star-shaped marks indicate rapid expansion events at MIS 5e (red), early MIS 3 (orange), MIS 1 (blue), and mid Holocene (green). c A geographic map showing an ancestral reconstruction of RrC, B. bengalensis, and M. musculus eastward dispersals over the last 130,000 years across the Indo-Malayan region ( Table 2). In RrC, the expansion event was likely linked to the dispersal event from Myanmar to Bali, Indonesia (Aplin et al. 2011) distribution of B. bengalensis may be explained as a recent invasion from another geographic area, namely the Indian subcontinent. This possibility may also explain the apparent geographic partitioning of haplotypes; haplotypes are spatially separated from one another (Mori et al. 2020). Contrary to the case of R. exulans, the dataset of B. bengalensis haplotypes from Myanmar yielded a rapid expansion signal, with an estimated expansion initiation time of 4000 (CI: 1900-6700) years ago. This implies that a geographic area including Myanmar experienced noteworthy historical events at this time (Fig. 4), and a plausible explanation is the prehistoric development of paddy rice cultivation in the geographic area (Fuller 2011). However, one haplotype from Pyay (Locality 4 in Fig. 1b) was markedly divergent from the other haplotypes, with genetic distance of 0.011-0.012 (Fig. 2c). This may imply that the species has been distributed in Myanmar for a period of time longer than 100,000 years. Alternatively, it could be the result of a separate introduction from a different location. The small number of samples analyzed and the lack of sequence data from the Indian subcontinent are issues that we plan to address in future research to clarify the evolutionary history of this species.
For M. musculus, all haplotypes from Myanmar belonged to CAS-1, which is known to occupy Pakistan, India, and Southeast and East Asian countries including China (Suzuki et al. 2013;Chingangbam et al. 2015;Bibi et al. 2017;Myat Myat Zaw et al. 2019). In this study, Myanmar haplotypes, together with those of the neighboring area of easternmost India (Manipur), are shown to form a star-like cluster ( Table 2, Fig. 3, Supplementary Fig. S1), implying that, for commensal mice, the Arakan Mountain system has no role as a geographic barrier between India and Myanmar. The appearance of a signal of rapid expansion may suggest that a single colonization event from either India or China, both of which are considered centers of diversification of the CAS-1 system (Suzuki et al. 2013;Kuwayama et al. 2017), shaped the genetic diversity of the CAS-1 mice of Myanmar. The expansion event for the CAS-1 sequences from Myanmar and its neighboring area (n = 19) was estimated to have occurred 5000 (CI: 500-7900) years ago (Table 2). These results suggest that human activities during the Holocene were involved in driving the sudden expansion of M. musculus in Myanmar (Fig. 4), as has been predicted previously (Myat Myat Zaw et al. 2019). In our previous study, we showed that the CAS-1 lineage experienced two expansion events (Suzuki et al. 2013;Kuwayama et al. 2017), which occurred 7000-8000 years ago in India and around 4000 years ago in China. Given the suggestion that rice fields emerged in eastern India and Yunnan, China, approximately 4000-4500 years ago (Fuller 2011), it is plausible that the sudden expansion event in Myanmar was triggered by the spatial expansion of paddy fields.

Conclusion
We assessed the evolutionary histories of murine species commonly seen in human-occupied areas of Myanmar, by comparing the sequence variations of Cytb in the Pacific rat (R. exulans), roof rat (RrC), lesser bandicoot rat (B. bengalensis) and house mouse (M. musculus). Cytb sequences of RrC, R. exulans, and M. musculus from Myanmar were determined and used for phylogeographic analysis, along with sequences obtained from databases. In roof rats, the Myanmar Cytb sequences were shown to belong to RrC Lineage II, which ranges from northeast India to Thailand, Vietnam, China, and Japan. A majority of the haplotypes from Myanmar and the neighboring region of Manipur, India can be integrated into two nested clusters containing signals of rapid population expansion. Employing a time-dependent evolutionary rate, the expansion start times based on nucleotide substitution levels were calculated to be 63,000 (CI: 44,000-77,000) and 140,000 (CI: 101,000-183,000) years ago, suggesting that both the colder period and subsequent warmer period during the transition from MIS 4 to MIS 3 and MIS 5e were involved, respectively (Fig. 4). This finding implies that species found in subtropical areas had been affected by late Quaternary environmental fluctuations, as had those in areas of higher latitude Hanazaki et al. 2017;Honda et al. 2019).
The mtDNA of RrC Lineage II exhibits a higher level of genetic diversity in Myanmar, implying that the area is its point of origin. The genetic diversities of R. exulans, M. musculus, and B. bengalensis from Myanmar, in contrast, were relatively low, suggesting that the three species were introduced to Myanmar from neighboring countries relatively recently. The haplotypes of B. bengalensis and M. musculus from Myanmar exhibited population expansion signals with an estimated expansion time of 4000-5000 years ago (Fig. 4), probably linked with prehistoric rice field development. By contrast, no expansion signal associated with mid-Holocene human activities was detected for R. exulans or RrC in Myanmar, suggesting that rice field development may not have been significantly associated with the population dynamics of these two rat species. In a future study, we will address the lineage dynamics of these commensal rodents to obtain a comprehensive view of their evolutionary histories at higher resolution using whole mitochondrial genome sequences (Li et al. 2020). Furthermore, we wish to acknowledge the two anonymous reviewers for their constructive comments. This study was supported by a grant-in-aid for Scientific Research (B) to HS from the Japanese Society for the Promotion of Sciences (JSPS; 17H04604). This study was conducted in part with the support of a JSPS KAKENHI Grant Numbers (JP18H05508, 24405045) and grants-in aid from the Japan Agency for Medical Research and Development (AMED) (JP15fk0108005, JP16fk0108117, JP17fk0108217, JP18fk0108017, JP19fk0108097, and JP20fk0108097).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.