Natural Diversity and Phylogeny of Asian Red-Cheeked Squirrels (Rodentia, Sciuridae, Dremomys) in Eastern Indochina

Based on new molecular data for mitochondrial (Cyt b) and nuclear (IRBP, RAG1) genes, as well as an extensive analysis of morphological material, we accessed actual species taxonomy and relationships among Asian red-cheeked squirrels Dremomys distributed in eastern Indochina and southern China. Phylogenetic analyses demonstrated that Asian red-cheeked squirrels, which are currently attributed to D. rufigenis, are not homogenic but instead consisted of two independent species-level clades—northern and south-central. The latter clade was additionally subdivided into two highly divergent clades based on Cyt b gene phylogeny. In spite of multidimensional statistics approach applied (PCA) only minor cranial differences were found between populations of study what lay a basis to treat it as cryptic species. Based on our findings, red-cheeked squirrels inhabit northern Vietnam and southern China, which are usually attributed to D. rufigenis, should be treated as distinct genetic species D. ornatus Thomas, 1914. In ones turn, based on its peculiar external morphology we can attribute the specimens from southern and central Vietnam to D. rufigenis proper and treat them as D. rufigenisfuscus Bonhote, 1907 and D. r. laomache Bonhote, 1921, respectively.

A number of different forms have been proposed as subspecies for the Asian red-cheeked squirrel D. rufigenis. Only three subspecies were recognized by Moore and Tate (1965): D. r. rufigenis Blanford, 1878, D. r. belfieldi Bonhote, 1908 andD. r. opimus Thomas et Wroughton, 1916, with two more taxa (fuscus Bonhote, 1907 and ornatus Thomas, 1914) distinguished by later reviews (Corbet and Hill, 1992;Thorington and Hoffmann, 2005); with opimus has been treated as junior synonym for adamsoni Thomas, 1914. According to recent review two subspecies are distributed in Vietnam: D. r. rufigenis inhabits northern Vietnam and Yunnan, and D. r. fuscus occurs in southern Vietnam and Laos (Thorington et al., 2012). The main object of this report was to investigate morphological and genetic diversity of D. rufigenis in Vietnam to delineate both taxonomical content and morphological variability of this taxon.

MATERIALS AND METHODS
A number of Dremomys specimens, has been obtained in Vietnam during a series of field expeditions of the Joint Russian-Vietnamese Tropical Research and Test Centre between 1998 and 2018. All these materials were obtained in full agreement with actual Vietnamese regulations in the field of Nature Protection and Biodiversity Conservation. We have examined 64 skulls, 107 skins, and 7 ethanol-preserved museums specimens of D. rufigenis, as well as 12 skulls and 13 skins of D. pernyi which constitute the largest data set ever investigated for this species within single survey. In general, specimens were obtained from 47 localities in Indochina and Southern China ( In total, 58 intact skulls of adult D. rufigenis from 20 localities in Vietnam (see Supplementary Table S1) were measured for morphometric comparison and analysis. Age was assessed using molar eruption and cranial seams conditions. For comparison, we used 9 skulls of D. pernyi obtained from the collections of ZIN, IEBR and ZMMU. Nineteen craniodental measurements were taken for each skull using digital calipers to the nearest 0.1 mm: greatest length of skull (ONL), palatal length (PL), interorbital breadth (IB), braincase breadth (BBC), braincase height (HBC), zygomatic breadth (ZB), interorbital breadth (IB), supraorbital breadth (SB), diastema length (LD), length of the nasal bones (NL), palatal length (LBP), postpalatal length (PPL), upper molar row length (CLM), lower molar row length (CLm), breadth across the palatal bridge at the level of the first molar (BBP1), breadth across the palatal bridge at the level of the third molar (BBP3), length of incisive foramina (LIF), breadth of the mesopterygoid fossa (BMF), A one-way analysis of variance (ANOVA) was performed to test the differences between groups for all cranial variables and indices. A principal component analysis (PCA) was computed to evaluate the diversity within Dremomys populations and species. The software package Statistica 8.0 (StatSoft Inc., Tulsa, OK, USA) was employed for all statistical procedures.
Total DNA from 96% ethanol-preserved muscle tissue was extracted using a routine phenol-chloroform-proteinase K protocol (Kocher et al., 1989;Sambrook et al., 1989). The DNA was further purified by twofold ethanol precipitation or a DNA Purification Kit (Fermentas, Thermo Fisher Scientific Inc., Pittsburgh, PA, USA). Four genes, both mitochondrial and nuclear, which proved to be useful for the phylogenetic analysis of Asiatic Sciuridae (Oshida et al., 2000(Oshida et al., , 2006Mercer and Roth, 2003;Steppan et al., 2004;Chang et al., 2010;Li et al., 2008Li et al., , 2013Hawkins et al., 2016) were targeted. Three of these genes were used in all analyses, including the complete cytochrome b (Cyt b) gene (1140 bp), a portion (up to 1610 bp) of the first exon of interphotoreceptor retinoid-binding protein (IRBP also known as RBP3), and a portion (1571 bp) of the first subunit of the recombination activation factor gene (RAG1). For some samples, we also sequenced the 5′-proximal 680bp portion of subunit I of the cytochrome oxidase gene (COI).
The Cyt b gene was amplified using the primers L14723 and H15915 (Irwin et al., 1991) and COI gene was run using the universal conservative primers BatL 5310 and R6036R (Kocher et al., 1989;Irwin et al., 1991). The following universal PCR protocol was used to amplify both mtDNA fragments: initial denaturation for 90 s at 95°C, denaturation for 30 s at 95°C, annealing for 60 s at 52°C, and elongation for 30 s at 72°C, followed by terminal elongation for 2 min at 72°C. The PCR reaction was performed in a volume of 30-50 μL containing 2.5 to 3 μL 10× standard PCR buffer (Fermentas, Thermo Fisher Scientific Inc., Pittsburgh, PA, USA), 50 mM of each dNTP, 2 mM MgCl 2 , 10-12 pmoles of each primer, 1 U of Taq DNA polymerase (Fermentas, Thermo Fisher Scientific Inc., Pittsburgh, PA, USA), and 0.5 μL (20-50 ng) of total DNA template per tube. The IRBP gene was amplified using the IRBP125f, IRBP223, IRBP1435r, IRBP1125r, and IRBP1801r primers (Suzuki et al., 2000) according to the method of Stanhope et al., (1992). A fragment of the RAG1 gene was obtained using the previously published primers S70 and S71; these primers were used both for the PCR and for sequencing as described in Steppan et al. (2004).
The PCR products were purified using a DNA Purification Kit. The double-stranded DNA products were directly sequenced in both directions using an ABI PRISM 3730xl Genetic Analyzer (Applied Biosystems, USA) and the BigDye® Terminator v3.1 Cycle Sequencing Kit in agreement with the manufacturer's protocol. All sequences were deposited in Gen-Bank (KX171231-KX171276).
In total, seventeen animals from eight geographical localities were genotyped: 17 for Cyt b, 10 for IRBP, 12 for RAG1, and 4 for COI (see Supplementary Table S2). For comparison, we analyzed all the Cyt b, IRBP and RAG1 gene sequences from the Dremomys species that are currently available in the GenBank/NCBI/JNDB databases and a few from the Callosciurus species used as an outgroup to root trees (Supplementary Table S2).
The sequences were aligned using BioEdit v 7.0 (Hall, 1999) and Clustal W (incorporated into BioEdit and MEGA 6) software and were verified manually. The basic sequence parameter calculations (i.e., variable sites, parsimony informative sites, base composition biases, nucleotide frequencies, and nucleotide substitution tables), the best-fitting gene evolution models and inter-and intrapopulation divergence evaluations were performed using the MEGA 6 software (Tamura et al., 2013). The most frequently used algorithms, such as minimum evolution (ME), maximum parsimony (MP), and maximum likelihood (ML), were applied for the phylogenetic inferences. The phylogenetic trees were constructed based on three genes: Cyt b (1-1140), RAG1 (1141-2712) and IRBP (2713-3824 nucleotide positions) and for concatenated dataset. Nucleotide evolution models were tested using MEGA 6 model's module. The HKY + G (Hasegawa-Kishino-Yano, gamma parameter =0.16) substitution model was used for the Cyt b gene, whereas model K2 (Kimura 2 parameter) is proved optimal for the RAG1 gene; the more complex GTR + G + I (General Time Reversible, gamma parameter =0.21) model was better for IRBP and was also used for concatenated dataset. The gamma shape parameters were evaluated and calculated directly from a general dataset. Because considerable saturation was revealed, COI sequences were not used directly for phylogenetic reconstructions and are only presented here as genetic vouchers for DNA barcoding and genetic distance comparison. The robustness of the trees was assessed using a bootstrap procedure with 1000 replicates. Trees were constructed and visualized directly using MEGA 6 or TreeView 1.6.6 software (Page, 1996). Genetic distance between individuals was calculated by pairwise distance based on Cyt b and COI sequences using the maximum composite likelihood models (Tamura et al., 2004) under the most complex T3P (Tamura 3 parameter) algorithm that was implemented in MEGA 6. Tajima's Relative Rate Test was applied to test the equality of evolutionary rate between evolution lineages as implemented in MEGA 6 and resulted in no reasons to reject the null hypothesis of equal rates between lineages. The dating was calculated based on reference points for some bifurcation events calculated by Hawkins et al. (2016), under ML Clock Test dialogue (Tajima, 1993) in MEGA 6 which was applied to Cyt b Maximum Likelihood tree.
To clarify the genealogical relationships of Dremomys, a minimal spanning network was constructed based on Cyt b gene haplotypes to investigate the evolution of the lineages. We estimated the evolutionary relationships by constructing a maximum-parsimony network of Cyt b haplotypes using the Median-Joining method (Bandelt et al., 1999) as implemented in Network 5.0 software (http://www.fluxus-engineering.com/sharenet.htm) with the default parameters (Polzin and Daneschmand, 2003) and without MPrelaxation of the network.

Phylogenetic analysis.
We found out high genetic divergence among Dremomys species and geographic populations in Indochina (Supplementary Table S3). For the Cyt b gene, the distances between haplogroups that were found in the northern, central and southern Vietnam (11.5-19.3%, T3P) fell within the range of interspecific distances observed in rodents (Bradley and Baker, 2001).
Phylogenetic trees for the Cyt b, RAG1, and IRBP genes are presented in Figs. 2a-2c. As can be observed on the trees, a complex phylogenetic structure appeared within Dremomys, and the most spectacular pattern shows the eastern Indochinese populations. Based on the Cyt b gene, there are three well-supported and deeply divergent phylogenetic lineages of D. rufigenis (Fig. 2a); the positions of other branches are not always reliable. We label one of main branches as Northern (N) and two others, which appeared as sisters to the first branch, as Southern (S) and Central (C). The support levels for three Vietnamese lineages appeared reliable using all phylogenetic methods applied. Nevertheless, lower levels of bootstrap values were noted for the Chinese clades attributed to D. per-  nyi, D. lokriah and D. pyrrhomerus, possibly indicating the basal radiation of these phyla.
The GenBank samples NC026442 and KC447304 (unknown localities) constitute the most unstable branch, and this lineage may appear as a sister to Vietnamese D. rufigenis and to Chinese D. pernyi, or even to D. pyrrhomerus, depending on the clustering method applied, but with low node support. Three Vietnamese branches of D. rufigenis corresponding to N-, C-and S-branches ( Fig. 2b) are supported by the IRBP gene with rather low reliability too. The structure of IRBP gene tree is more complicated, apparently because this gene is the most slowly evolve and the least variable compared with other genes used for mammalian phylogeny (Stanhope et al., 1992). The RAG1 gene tree is the most intriguing, as is illustrated in Fig. 2c. Despite the insufficient level of reliability observed (which most likely occurred because the sequences available for Chinese populations of D. pernyi were very short), the tree topology remained robust irrespectively to the methods, calculation algorithms and basic parameters applied. The northern Vietnamese populations clustered together with D. pernyi, whereas other populations constitute independent sister groups.
Faced with the situation where deeply divergent phylogenetic lineages were evident but their relationships could not be reliably supported by single-gene phylogeny, we performed an analysis based on concatenated dataset to improve fidelity and reliability of resulting tree. Combined analyses based on the three genes resulted in a well-supported tree (Fig. 3), where two major phylogenetic clades (N group and united S + C group) were evident for D. rufigenis in Vietnam.
The analysis of the median-joining network structure also shows that Vietnamese Dremomys consists of two independent subpopulations (Fig. 4). The first subpopulations, exactly N group samples constitute the torso of the network and tend to be closest to the Chinese D. pyrrhomerus. The central and southern Vietnamese populations constitute another lineage, which can be reliably connected to the torso only through the Southern haplotypes; the Central group haplotypes are much more distant. This pattern is inconsistent with the current species distribution pattern, but inevitably indicates that the southern and central Vietnamese populations should be regarded as independent entities from the Northern group.
Eastern Indochinese lineages appeared as separated from united generic torso, demonstrating an independent evolution process. The evolution trajectory of northern clade exhibits a clear longitudinal pattern and long branches. An evolution of S-C clade appears to be a progressive process from the southern to more northern haplotypes of this group, which can illustrate its natural history during the spread from south to north.
The network also clarifies the position of some northern branches of D. rufigenis and D. pyrrhomerus, showing that the phylogenetic tree is unsustainable  ( Fig. 2A). D. rufigenis appeared as the most divergent clade originating from the general pull of the Northern phylogenetic lineage. The considerable branch length of this clade is consistent with its specific status, whereas some other shorter branches can be regarded in context of intraspecific diversity.
Analysis of external characters. Pelage coloration and some other fur characteristics are subjected to considerable individual variation in Southeast Asian red-cheeked squirrels (Dao, 1969;Cao, 1984;Dao and Cao, 1990). This is most apparent for the populations from Thai Nguyen, Vinh Phuc, Nghe An, Quang Ninh and Lang Son provinces. Ventral side of animals collected in the same geographical area may vary from bright white to dirty-yellow or grayish color. Dorsal side and cheek colorations are also very variable. In some specimens the head, back and hips are lightbrown colored with buffy hue, whereas other animals have prominent brownish pelage, with reddish hue. These patterns may well be depending on seasonal variation, as most of light-colored specimens were collected from May to July.

Animals from central Vietnam (Thanh Hoa, Ninh
Binh and Nghe An provinces) and central Laos are appreciably lighter colored, with buffy-reddish and brownish hues of dorsal pelage. Their coloration is in concordance with that of laomache, a peculiar pallid race described from Laos (Thomas, 1921).
The light patches behind the ears and the pelt texture are the most distinctive features of various D. rufigenis populations. The populations from the N group have yellow or ochreous ear patches, which could be seen very well from the lateral view. The animals from southern and central populations have whitish ear patches which could be only seen clearly when looking from behind or above. It should be noticed that these white patches behind the ears were emphasized in the original species description for D. rufigenis (Blanford, 1878) and are also listed for Malayan (Bonhote, 1908) and Burmese (Thomas and Wroughton, 1916;Moore and Tate, 1965) specimens. D. fuscus from Nha Trang, Southern Vietnam has the same coloration pattern (Bonhote, 1907). On the other hand, buffy patches behind the ear were noticed for the Chinese D. ornatus and D. adamsoni from the northern Burma by Thomas (1914).

Analysis of cranial characters.
Means, ranges and standard deviations of the craniodental measurements and indices for three groups that were identified by the genetic analyses are given in Supplementary Table S4. As revealed by ANOVA, specimens of the Southern group significantly differ from the Northern and Central groups; many cranial characters (ZB, BBC, LB, HB, LB, HB, CLM, PL, and CLm) were smaller, and some indices (IZB and ILB) also significantly differ. In contrast, differences between the Northern and Central groups were low; the only significant difference among these groups was found in the breadth  The complex morphological variation in 58 Dremomys skulls was analyzed using a PCA based on all 30 craniodental measurements and indices. Most characters were highly positively correlated with the first component PC 1 (Supplementary Table S5), suggesting that PC 1 reflects differences in general size of cranium. The second component PC 2 was strongly correlated with the length of the auditory bulla (LB). The specimens of D. pernyi and D. rufigenis s. lato form two essentially discrete groups (Fig. 5). However, even multivariate analysis did not find any visually perceptible cranial characters suitable for direct identification of studied populations of D. rufigenis s. lato without special statistical procedures.

DISCUSSION
Taxonomical implications. Skulls of D. rufigenis from a few Indochinese localities were examined by Endo et al. (2003). They stated that the squirrels in the southern area (Malaysia) were found considerably larger than those in the northern areas (Southern Vietnam, Laos and Thailand; for example, maximum skull length 51.47 mm, SD = 1.24 against 50.32 mm, SD = 1.44 and 48.9 mm, SD = 2.32, respectively for adult males). The splanchnocranium was longer in the Malayan population and the interorbital space was narrower for the Vietnam-Laos and Thailand popula-tions. However, lacking any genetic data available, no taxonomical conclusions have been made based on these findings.
Based on our results, we must conclude that squirrel inhabiting northern Vietnam and southern China, which is usually regarded as D. rufigenis, can be treated as distinct cryptic species D. ornatus Thomas, 1914, rather than a subspecies of D. rufigenis. As was mentioned above, molecular diversity among studied populations was found much higher than morphological one. The situation is not so surprising. There are many examples of pairs or groups of species that, while may demonstrate lack of morphological differences, remain valid species. These are cryptic or genetic species, which have been recognized recently mainly or even solely based on drastic genetic segregation (Baker and Bradley, 2006). Some mammalian species-Loxodonta cyclotis (Roca et al., 2001), Carollia sowelli (Baker et al., 2002;Solari and Baker, 2006), Neotoma macrotis (Matocq, 2002), Myotis occultus (Piaggio et al., 2002), Notiosorex cockrumi (Baker et al., 2003), Lophostoma equatorialis (Baker et al., 2004), Peromyscus schmidlyi (Bradley et al., 2004a), Reithrodontomys bakeri (Bradley et al., 2004b), Lagidium boxi (Spotorno et al., 2004) and Thylamys cinderella (Braun et al., 2005)-have been recognized, described or justified using DNA sequence data alone. Many other examples described based on other genetic markers are known within Sicista, Microtus, Sylvaemus (see Sokolov et al., 1981Sokolov et al., , 1986Sokolov and Baskevich, Fig. 5. Ungrouped morphometric separation (principal components analysis) of Dremomys specimens from Eastern Indochina; the data were drawn from craniodental measurements and indices (see Supplementary Table S4). Symbols: pern-D. pernyi, N, C, and S-specimens of D. rufigenis s. lato from the Northern, Central and Southern groups, respectively. 1988; Chelomina et al., 1998;Conroy and Cook, 2000;Hille et al., 2002).
In the situation we faced, it is necessary to provide here formal description to reestablished D. ornatus as independent taxon in agreement with its initial description of Thomas (1914). Thomas, 1914  Description. According to Thomas (1914: 26): "General appearance of rufigenis, but skull characters as in belfieldi. Coloration almost as in true rufigenis, but the back is clearer olive without any muddy or drabby tinge, and the rump with scarcely any sufficient of brown. Undersurface as rufigenis. Sides of face brighter ferruginous, the rusty extending further back on the upper surface of the muzzle, nearly to the level of the ears. Patches behind ears buffy. Skull very different of that of rufigenis, but indistinguishable from that of belfieldi. Larger then in rufigenis, the muzzle very long pinched in at sides; the nasals long, their hinder ends directed transverse. Forehead broad, flat with more transversely directed posterior processes. Brain case large. Zygomata widely expanded. Hindfoot, 44. Skull: greatest length 58.2; condylo-incisive length 50.3; zygomatic breadth 31.6; nasals 21.7; interorbital breadth 16.7; breadth of brain case 25; palatial length 24.2; p 4 and three molars 9.7." Diagnosis. "Distinguished from rufigenis by its large skull and long muzzle and from belfieldi by its brighter coloration" (Thomas, 1914). The species is well differentiated from other Dremomys based on mitochondrial (Cyt b, COI) or nuclear (RAG1, IRBP) genes. The divergence level from other species-level lineages is d T3P 0.187-0.269 for Cyt b and 0.0418-0.047 for COI genes. The DNA sequences deposited in GenBank (KX171244-KX171247; KX171258-KX171260; KX171264 and KX171273-KX171276) can be regarded as genetic vouchers for the species.

Dremomys ornatus
Comparison. D. ornatus is cryptic genetic species, sibling with D. rufigenis. General coloration pattern is rather variable and depends on the season and local variation. The most prominent external characters are the color of the spots behind the ears and the pelt texture. D. rufigenis have the ears scarcely covered by hairs; the patches behind the ears are brightly white and can be only seen from behind or above (Fig. 6). In contrast, D. ornatus has well-haired ears, the ochreous ear patches can be well seen from the lateral side (Fig. 7). The skull (Supplementary Figs. S1-S2) is usually larger in D. ornatus (ONL up to 58 mm or more), but no distinct morphological cranial characters currently traced to distinguish it reliably from D. rufigenis by skull morphology.
Distribution. Genetically confirmed specimens of D. ornatus found in northern Vietnam (Ha Giang, Cao Bang and Son La provinces). Based on the morphological data, the species distributed throughout northern Vietnam, northward to Yunnan, and probably, to Guangxi provinces of China. D. ornatus meets the congeneric species, D. rufigenis laomache at Son La and Yen Bai provinces; thus, the southern limits of its distribution lie between 17°30′ and 17°45′ N. There is no genetic evidence of sympatry, but the presence of unusually colored specimens at the area of possible species contacts leaves room to suspect natural hybrid-  ization. Tail of these animals is very light colored, from ochreous-yellow to whitish, on its ventral side (Fig. 8), whereas it is orange-red in most of populations of this species. Such unusually colored individuals, who appear at notable part of populations of red-cheeked squirrels in central Vietnam, have even prompted some authors (Dao, 1969) to attribute these individuals to Perny's squirrel, D. pernyi flavior, what is an obvious mistake.

Notes on D. rufigenis sensu stricto composition.
Judging by the level of genetic diversity, both C-and S-lineages should be attributed to D. rufigenis, these groups can be treated as distinct subspecies. In agreement with ICZN (1999), there are only two appropriate nomena available, exactly laomache and fuscus (see Fig. 1, where all type localities for Indochinese redcheeked squirrels are indicated).
As stressed by Thomas (1921: 182-183), pale color is the main characteristic feature of D. r. laomache. Thomas wrote as follows: "A pallid strongly grizzled race of D. rufigenis. General characters about as in true rufigenis, the skull of similar proportions not excessively elongated, the postauricular patches not especially conspicuous, the under surface of about the same degree of greyness. However, the dorsal color is paler and greyer and much more coarsely speckled, the buffy subterminal lings on the hairs over 2 mm in length compared with rather less than 1 mm in rufigenis and its other subspecies. Rufous of face, rather brighter than in rufigenis. Postauricular patch not large or very conspicuous, its hairs white with their ends buffy, more broadly buffy than in rufigenis. … [This species] readily distinguishable by its pale and comparatively coarsely speckled back from any of the races hitherto described. D. r. fuscus, its near neighbor in Annam, is a dark-colored form with conspicuous white postauricular patches." According to original description, no morphological features other than coloration pattern were reported (Thomas, 1921), that is in full concordance with our result of statistical analyses. Of note, the descriptions presented are in complete agreement with the coloration pattern of the Laotian specimen (ZIN 98317), which collected near the type locality of D. r. laomache. Most of central Vietnamese specimens cannot be called as completely pale, although many of them are truly appreciably lighter then both southern D. r. fuscus and northern D. ornatus. Thus, the presence of lightcolored individuals is characteristic of the C-lineage, which we suppose to attribute to D. rufigenis laomache Thomas, 1921. At ones turn, S-lineage can be associated with D. r. fuscus, whose subspecific rank is generally accepted and genetic attribution is obvious.
Undoubtedly the taxonomic content of D. rufigenis is required the further revision. Lacking of appropriate sampling, it is currently difficult to appreciate the limits and range of nominotypical D. r. rufigenis regarding its relationship with the morphologically distinct southern Malayan populations (Endo et al., 2003) and northern D. lokriah. Perhaps, the former populations   can be attributed to D. r. belfieldi Bonhote, 1908 described from Selangor, Malaysia. Distribution and natural area formation. The redcheeked squirrels are strictly confined to forest environments, and their dispersal is restricted by the timberline. They cannot move beyond these limits or cross wide deforested areas. The main natural event that likely contributed to the dispersal, population segregation and speciation of squirrels is, most probably, repeated disjunction-reconnection events of natural populations that were associated with oscillation of areas covered by tropical forest during the late Miocene-Holocene (Hall, 1998).
To date, the most ancient fossils representative of the genus, D. primitivus Qiu, 2002, known from Lufeng, Yunnan, and have been dated as recent as 8-9 Mya (Qiu, 2002). According to Hawkins et al. (2016), the most recent common ancestor for all longnosed squirrels was estimated to live at 5.97 Mya (95% HPD 4.34-7.95), and the most recent speciation event was dated to 2.9 Mya (95% HPD 1.77-4.17, between D. gularis and one of the Chinese lineages of D. rufigenis s. lato); the average age of all recent species estimated as 4.4 Mya (range 5.97-2.9 Mya). Based on these evaluations the divergence event between the original south-central Vietnamese and northern Vietnamese populations occurred at 2.41 Mya (range 1.75-3.08 Mya). The split time between the Southern and Central phylogroups was evaluated to occur 1.80 Mya (range 1.3-2.31 Mya). This timing coincides with two most recent and greatest sea-level regressions (up to 60-70 m below recent sea level) to occur during the Pleistocene (Haq et al., 1987). This fact support the longitudinal isolation and suspect the specific status of these populations.
During the Pleistocene glaciations, the forest edge was located at lower elevations (Cox and Moore, 2000), and substantial forest contraction was discovered to exist during the last glacial periods, as evidenced in Peninsular Malaysia and Palawan (Wurster et al., 2010). Cranbrook (2000) and Gathorne-Hardy et al. (2002) stressed that most areas in this region were covered by savannahs during long periods of the Quaternary.
Recent studies on the biogeography and paleoenvironments of the Sunda Shelf have suggested that the interactions between climate and sea level and their effects on the distribution of the fauna and flora is much more complex (van den Bergh et al., 2001;Meijaard and Grove, 2006). Many islands within the Sunda Shelf were connected to each other by forested environments (Anshari et al., 2000;Kershow et al., 2001;Meijaard, 2003;Cannon et al., 2009;Woodruff, 2010;Raes et al., 2014;Leonard et al., 2015). Similarly, continuous forested cover in Indochina existed during a considerable period of the Pliocene and Pleistocene (Meijaard and Groves, 2006), enabling the direct contact, dispersion and genetic exchange of western and eastern Indochinese populations, which are currently interrupted by huge deforested areas in central Thailand and western Cambodia. No comparative genetic data regarding Malayan populations are available to sustain this proposal; nevertheless, based on data about tropical forest persistence and the distribution of forest-dependent species on islands of the South China Sea and another arboreal species distribution pattern (Meijaard et al., 2003;Balakirev et al., 2021), the connection between southern Vietnam and Malacca may even be more probable than that that with Peninsular Thailand. The complex paleogeographic history of the region resulted in the remarkable biodiversity and complex intraspecific structure of many species observed in Indochina, southern China and Sundaland (de Bruyn et al., 2014;Abramov et al., 2014;Leonard et al., 2015;Hu et al., 2021). During hot and arid periods, forested areas persisted only in the highlands. Isolated populations of Dremomys apparently inhabited these forest refugia during very long historical periods. Environmental barriers, such as the above-mentioned savannah-like areas, repeatedly appeared in the lowlands and reduced the genetic exchange between populations inhabiting isolated highland mountain forests, thus favoring genetic diversification and speciation events. The climate of this period sustained the persistence of moist mountain forests in Dalat, Kon Tum and more northern localities along the Annamite Range, such as the Bolaven, Phouane, and Nakai plateaus. The orography of Yunnan is even more complex and favors the speciation process. Perhaps, further studies on the genetic diversity across the Dremomys will uncover further undescribed squirrels' species in the mountains of Southern China and the Malay Peninsula, as has been found both for this genus (Li et al., 2008) and for related groups such as Tamiops (Chang et al., 2011) and Petaurista (Li et al., 2013).
The tectonic movements and alteration between glacial and interglacial periods during the Pliocene, along with the further uplift of the Himalayas approximately 2-3 Mya (Zheng et al., 2000), might have been the major factors responsible for the dramatic genetic variability and sympatric distributions of many Dremomys species, such as D. ornatus, D. gularis and D. pyrrhomerus, in southern China. Recent populations of D. r. fuscus and D. r. laomache might originated due to such favorable interactions. As evidenced by genetic and morphological data, these populations may originally descend from western Indochinese, not Chinese populations, which they only met much later due to secondary contact during the recent natural range renewal. Similar events were assumed earlier for some other arboreal species (Balakirev et al., 2012(Balakirev et al., , 2021Meschersky et al., 2016;Balakirev and Rozhnov, 2019).