A Remarkable Case of Micro-Endemism in Laonastes aenigmamus (Diatomyidae, Rodentia) Revealed by Nuclear and Mitochondrial DNA Sequence Data

L. aenigmamus is endemic to the limestone formations of the Khammuan Province (Lao PDR), and is strongly specialized ecologically. From the survey of 137 individuals collected from 38 localities, we studied the phylogeography of this species using one mitochondrial (Cyt b) and two nuclear genes (BFIBR and GHR). Cyt b analyses reveal a strong mtDNA phylogeographical structure: 8 major geographical clades differing by 5–14% sequence divergence were identified, most of them corresponding to distinct karst areas. Nuclear markers display congruent results but with a less genetic structuring. Together, the data strongly suggest an inland insular model for Laonastes population structure. With 8 to 16 evolutionary significant units in a small area (about 200×50 km) this represents an exceptional example of micro-endemism. Our results suggest that L. aenigmamus may represent a complex of species and/or sub-species. The common ancestor of all Laonastes may have been widely distributed within the limestone formations of the Khammuan Province at the end of Miocene/beginning of the Pliocene. Parallel events of karst fragmentation and population isolation would have occurred during the Pleistocene or/and the end of the Pliocene. The limited gene flow detected between populations from different karst blocks restrains the likelihood of survival of Laonastes. This work increases the necessity for a strict protection of this rare animal and its habitat and provides exclusive information, essential to the organization of its protection.


Introduction
Eroded limestone outcrops form a prominent part of the landscape in Southeast Asia. Because of high limestone weathering [1,2] erosion in these areas resulted in scattered, isolated limestone hills with steep flanks called karsts towers. Despite a high diversity of habitat specialists and endemic taxa [3][4][5][6][7][8][9], these limestone karsts remain among the least studied ecosystems in Southeast Asia. Between 1985 and 2004 they contributed only to 1% of the global and regional biodiversity research output from terrestrial and freshwater ecosystems, while they cover around 10% of the land area in Southeast Asia [4]. However, in the Khammuan Province of Lao PDR, several new endemic vertebrate species were recently described: a bird, Pycnonotus hualon [3], a bat, Hipposideros scutinares [10], a gymnure, Hylomys megalotis [11], and a murid rodent, Saxatilomys paulinae [12].
Laonastes aenigmamus, recently described from this region by Jenkins et al. [13], is of special interest. Using comparative morphological and molecular data, this new species was initially described as a member of a new genus and ranked within the Hystricognathi rodents in a new family: the Laonastidae [13]. But, further studies demonstrated: first, that L. aenigmamus is a member of the Diatomyidae, a fossil family known from early Oligocene to late Miocene in Pakistan, India, Thailand, China, and Japan [14]; second, that the Diatomyidae are the sister group of the Ctenodactylidae, a family of small rodents found in rocky deserts across the northern parts of Africa; third, that together with the Hystricognathi, the Diatomyidae and the Ctenodactylidae form the suborder Ctenohystrica [15].
In a previous study we presented a preliminary analysis of the genetic diversity of L. aenigmamus [16]: 52 specimens were sampled and the population structure was surveyed by sequencing 887 base pairs of the Cytochrome b (Cyt b) gene. Phylogenetic and haplotypic network reconstructions revealed three well-supported and rather divergent lineages, suggesting that L. aenigmamus populations are geographically structured. However, altogether these samples represented a limited part of the estimated range of the species; most of them were sampled in local markets and their exact geographic origin was impossible to determine; this study was based on a fragment of a single mitochondrial gene.
Before its scientific discovery, the kanyou was considered as a kind of game, trapped and eaten by the villagers. Since 2008, the conservation of L. aegnimamus is regulated in Lao PDR and since 2009 this species is listed as «Endangered» on the IUCN Red List. During the last four years we add the unique opportunity to accompany local officers in the field when they explained this change to local populations, and to sample for molecular analyses the last specimens that were captured by traditional hunters. We were able to sample 137 individuals for which the exact locality of collect was known. It represents 38 localities probably spanning the whole geographical range of the species. Thus, the sampling obtained for this study is exceptional both in terms of number of individuals and number of localities sampled and will not be able to be carried out again. In the present study we: (1) improve evaluation of the genetic diversity and delineation of the phylogeographical structure of L. aenigmamus, by sequencing both nuclear and mitochondrial markers of a larger number of geographically referenced individuals; (2) obtain a better definition of the geographical distribution of the species; (3) highlight the spatial and temporal history of the limestone formation of the Khammuan Province, using the evolutionary history of this species as a guide. Our analysis strongly underlines the necessity for a strict protection of this rare animal and its habitat and provides exclusive information, essential to the organization of its protection. All the specimens used in this study were captured by traditional hunters. When an animal was captured alive all efforts were made to minimize suffering (animals were maintained quietly in appropriate caging with enough food and moisture and adequate temperature), and it was handled under the guidelines of the American Society of Mammalogists [17]. Alive animals were euthanized by the injection of a lethal dose of isofluorane, followed by cervical dislocation.

Specimens examined and mitochondrial DNA sequencing
Specimens whom geographical provenance was not precisely established were excluded. The present study includes 137 individuals representing 38 localities ( Fig. 1; Table 1; see Table  S1 for the list of specimens used in this study, Genbank numbers and Museum numbers). DNA was extracted from ethanolpreserved tissues by the CTAB method [18]. The Cyt b gene was amplified for 135 specimens using polymerase chain reaction (PCR) primers L14723 and H15915 [19]. PCR conditions were the same as those described in Nicolas et al. [20]. We also amplified and sequenced two nuclear gene fragments: the Intron 7 of the beta-fibrinogen (BFIBR) and the exon 10 of the growth hormone receptor (GHR). Amplification and sequencing of these two nuclear genes were performed for 83 specimens and 80 specimens for BFIBR and GHR genes respectively (for technical reasons amplification of some individuals for one gene or the other failed), representing all known mtDNA groups. For the BFIBR we used the primers BFIBR1 and BFIBR2 published by Seddon et al. [21], and an annealing temperature of 50uC. For the GHR we used the primers GHR1 and GHR2 published by Adkins et al. [22] and the primers GHR7 and GHR8 published by Lecompte et al. [23], and an annealing temperature of 55-60uC. The doublestranded PCR products were purified and sequenced at the Genoscope (Ivry/Seine, France). Sequences were aligned by eye. Sequences were sublitted to the Genbank database (HQ687370 to HQ687474, and JN796953 to JN797145).
The ''numts'' are usually characterized by the presence of indels, stop codons, frame-shift mutations and amplification of heterozygotes [24,25]. We have not observed any of these indications in our dataset. Moreover the base composition bias per codon position was not significantly different between individuals. The strongest compositional bias for all individuals was observed at the 3 rd codon position (A = 36%, C = 46%, G = 4%, T = 14%), and the least at the first codon position (A = 27%, C = 26%, G = 22%, T = 25%). This is congruent with the results obtained for other vertebrates species [26]. So we believe that numts were not present in our dataset
The MP analysis was performed with tree-bisection-reconnection (TBR) branch swapping option and 10 random addition replicates. We estimated the robustness of internal nodes by 1,000 bootstrapping replicates (each with a single replication of random addition of taxa). An equal weighting of character-state transformations was applied.
MrModeltest 3.04 [30] was used to evaluate the fit of 24 nested models of nucleotide substitution to the data. According to Akaike information criterion, MrModeltest recommended the GTR+I+G model for the Cytb dataset, GTR+G for the BFIBR dataset and HKY+G for the GHR dataset. These models were used in ML and Bayesian analyses. ML analyses were performed on the PHYML online web server [31] using the BIONJ distance-based tree as starting tree. Bootstrap analysis (500 replicates) was used to estimate the robustness of internal nodes.
In all MCMC analyses three heated chains and one single cold chain were employed, and runs were initiated with random trees. Two independent MCMC runs were conducted with six million generations per run; trees (and parameters) sampled every 100 generations. Stationarity was assessed by examining the average standard deviation of split frequencies and the Potential Scale Reduction Factor [32]. For each run, the first 25% of sampled trees were discarded as burn-in.
Following Huchon et al. [15] conclusions about Laonastes phylogeny, our Cytb and GHR phylogenetic trees were rooted with two members of the Ctenodactylidae family (Massoutiera mzabi and Ctenodactylus vali) and two Hystricognathi (Heterocephalus glaber and Thryonomys swinderianus). As no BFIBR sequences were available for these species, we rooted our BFIBR tree with another member of the Ctenodactylidae family (Ctenodactylus gundi) and one Hystricognathi (Heliophobius argenteocinereus).

Population genetics analyses
Haplotype diversity for each gene was represented using a haplotype network. Prior to this analysis, the existence of heterozygous positions for the two nuclear genes fragments was investigated and an input file constructed from this information using SeqPHASE [33]. The phase of each haplotype and its reconstruction were carried out using PHASE v.2.1.1 [34,35] by running the formerly built input file and by considering the default parameters of the software. This data file was used for the construction of a median-joining haplotype network, using the software NETWORK v4.500 [36].
We applied the SAMOVA procedure (Spatial Analysis of MOlecular VAriance, [37] to identify groups of geographically adjacent populations that were maximally differentiated based on sequence data. Using the computer program SAMOVA 1.0, we performed analyses based on 100 simulated annealing steps and examined maximum indicators of differentiation (F CT values) when the program was instructed to identify K = 2 through K = 14 partitions of populations. We analysed population structure with analysis of molecular variance (AMOVA) in Arlequin 3.5.1 [38]. AMOVA divides the total variance into additive components, that is, variation within populations, among populations within groups and among groups. A population was defined as all individuals coming from one geographical locality; and groups of populations were defined according to the results of the SAMOVA analysis.
In order to determine if each obtained mtDNA clade can be treated as a different population, we applied the Snn test [39] on each nuclear gene fragment to test for genetic differentiation between the obtained mt DNA clades. Calculation of the P-value of the Snn was performed with a permutation test with 1000 replicates implemented in the software DnaSP v.5 [40].
Where applicable, we estimated several statistics to describe and compare the major clades supported by our phylogenetic analyses. Number of haplotypes, number of polymorphic sites, nucleotide diversity, haplotype diversity and average number of nucleotide differences were calculated for all clades with more than 10 individuals using DNASP v5 [40]. We also applied the Tajima's D [41], Fu and Li F* and Fu and Li D* [42] neutrality tests to each nuclear and mtDNA lineage using the same software. The significance of these tests was calculated using 10,000 coalescent simulations.
The plausibility of an isolation-by-distance scenario was explicitly tested by performing Mantel's tests [43] using Arlequin 3.5.1 [38] to analyze the relationships between genetic (mean number of pairwise nucleotide differences) and geographical distances between sampling localities. The geographical distances between localities were calculated from GPS coordinates in the software GENALEX v6 [44]. Data were permuted 1,000 times to estimate the 95% upper tail probability of the matrix correlation coefficients.

Time of divergence
Using the clades identified in the phylogenetic analyses, a Bayesian MCMC approach was used to calculate the time since genes were last shared (the most recent common ancestor, TMRCA). TMRCA can be used as a proxy for ancestral population age. Clades were dated using BEAST v1.4.6 [45]. As a calibration point we used the split between the Ctenodactylidae and Laonastes, estimated at 44.363.5 Mya [15]. A Bayesian skyline coalescent model of population size was specified [46]. Divergence times and their credibility intervals were estimated using a relaxed clock model with branch rates drawn from an uncorrelated lognormal distribution. Two independent runs of 10 million generations each, with burning of 1 million generations, were performed and later, they were combined in Tracer version 1.4 [47], which provides options for examination of the effective sample size values and frequency plots, in order to check if the mixing of the MCMC chain was adequate. Samples from both runs were combined using the software LogCombiner v.1.4.7 (available in the beast package), and a consensus chronogram with node height distribution was generated and visualized using

mtDNA
A fragment of 983 bp of the Cyt b gene was retained for our final analyses. Among the 135 sequences analyzed, 55 different haplotypes with a total of 286 polymorphic sites can be identified. The average number of nucleotide differences is 50.56. The haplotype and nucleotide diversities are high: 0.95160.009 and 0.05160.002, respectively. The results of the MP, ML, Bayesian and network analyses are congruent ( Fig. 2 and 3) and highlight the existence of several well-supported clades within L. aenigmamus. Our phylogenetic analyses reveal a strong phylogeographical structure: specimens from a single clade always originate from adjacent localities on the map (Fig. 4) Phylogenetic resolution between clades is low. Only two wellsupported nodes are recovered, the first one grouping A1-A2-A3-A4-A5-B1-B2-C-D1-D2-E1-E2-E3-F and the second one grouping E1-E2-E3-F.
Average percentages of pairwise differences (Kimura 2-parameter) between clades varied from 5% between clades B and C to almost 14% between clades A and H ( Table 2). Average percentages of pairwise differences between subclades vary from 1% to 4.5% (Table 3). Within clades, average percentages of pairwise differences vary from 0 to 3%, and within subclades it varies from 0.1 to 1.8%.
Nucleotide and haplotype diversities of the main clades are given in Table 4. Neutrality tests were never significant, except the Tajima's D test for the D2 clade (Table 4).
Results from the SAMOVA analyses showed, as expected [37], that when the number of groups of populations (K) increased, the value of F CT increased. F CT only increases in small increments for K = 6 to K = 14, suggesting that adding extra groups only moderately improved the model of population structure. The group membership at K = 6 is highly congruent with the groups found in the phylogenetic and network analyses: these 6 groups correspond to clades A, B+C, D, E, F, G+H. Genetic structure was highly structured among these groups of populations (AMOVA; permutation tests highly significant, P,0.00001), meaning that most of the genetic variability was partitioned among groups (77.28%), and only a small proportion of the variability was partitioned among populations within groups (14.71%) or within populations (8.003%). The configuration defined for K = 6 is preserved for K = 7 to K = 14, meaning that barriers are incrementally added one after the other. Clades B and C are separated for K = 8, and clades G and H are separated for K = 9.
A significant correlation between geographical and genetic distance is observed (Mantel tests, P,0.0001, r = 0.625; Fig. 5) and, as expected, the genetic distance between geographically distant populations is always high. Conversely, high genetic differentiation variability may also be observed between neighboring populations. Figure 5 shows that for a small geographical distance (less than 39 km), the genetic distance can be either low (0 to 40; corresponding to intra-clade differences) or high (up to 118 pairwise nucleotide differences). Specimens coming from localities distant by less than 5 km from one another may differ by up to 58 nucleotide differences: between Phon Savanh North and Phon Savanh-km1. Whatever the geographical distance (19 to 153 km), specimens coming from the localities of Ban Vong Khon (clade G) and Ban Thalak (clade H) always differ by, at least, 90 nucleotide from all other specimens.
The TMRCA of all Laonastes sequences is estimated at 8.58 (range 5.84-11.61) Mya (Fig. 6). The TMRCA of clade E is 1.94 (1.20-2.79) Mya, and those of E-F is 3.20 (2.00-4.50) Mya. The TMRCA of most clades is recent (less than 1.5 Mya). nDNA A total of 83 specimens, including representatives from all the defined mtDNA lineages, were sequenced for the BFIBR gene. A total of 80 specimens, including representatives from all but one (no specimen of clade H could be amplified) of all the mtDNA clades were sequenced for the GHR gene. The obtained fragments have a length of 694 bp (34 polymorphic sites for the ingroup sequences) for BFIBR and 856 bp (41 polymorphic sites for the ingroup sequences) for GHR. Heterozygous specimens are found for both gene fragments. Overall nucleotide variation in the nuclear markers is low. This low level of sequence divergence results in unresolved phylogenetic trees ( Figure S1): Only clade H is supported for BFIBR sequences and clade G is supported for GHR sequences.
In the median joining networks clade H is well separated for BFIBR gene, and clade G is well separated for the GHR gene (Fig. 7). Moreover clades A1, A2, A, B1, C, D, D1, D2, E1, E3 and F all have one or several privative haplotypes for BFIBR, and the same is true for clades A1, A, B2, B1, D1, D2, D, E1, E2 and E3 for the GHR. When the two nuclear genes are combined (Fig. 7) resolution is better (only specimens for whom both genes were sequenced are included in the network, thus clade H is not represented): all mtDNA clades harbor one or several privative haplotypes (except clade A5), and clades A1-A2-A3-A4, D1-D2 and G are well separated.
The results of the Snn test (Table 5) indicate that the clades defined by mtDNA are significantly differentiated for both nuclear genes (except clades B and E). Subclade differentiation could only be tested between D1 and D2 and is significant for BFIBR only. All clades were therefore treated as different populations and demographic analyses were performed independently for those with appropriate sample size.
Nucleotide and haplotype diversities are given in Table 4. Neutrality test are never significant, except for the D1 clade and BFIBR.

Comparison of mtDNA and nDNA results
The high genetic divergence observed between mitochondrial lineages, is not fully corroborated by nDNA genes. However, given that: (1) private haplotypes are found for all clades, (2) Snn tests indicate that mtDNA defined clades are also significantly differentiated for both nuclear genes, (3) pattern of differentiation between clades is clearer when both nDNA genes are combined, and (4) between clade divergence events mostly occurred recently (during the Pleistocene); incomplete lineage sorting of the nuclear markers could explain the low genetic differentiation observed in nDNA. The higher evolutionary rate of mtDNA compared to nDNA, and its effective population size being one-quarter that of nuclear markers, allow a chance of recovering the pattern and tempo of recent historical events [49].
An alternative possibility would be that observed nuclear marker patterns could be due to a male-biased gene flow, with recurrent dispersion of male and rare dispersion of females leading to a contrasting pattern of nucleotide diversity between mtDNA and nDNA. Nothing is known on the mating and dispersion behavior in this species. However observations in the field, as well as metabolic and anatomical peculiarities of the species suggest that it is restricted to the karst habitat, with a low probability of survival if attempting to cross extensive open areas [50]. Thus male dispersion from one karst tower to another is doubtful, and it seems unlikely that a sex-biased gene flow could alone explain the contrasting pattern of differentiation observed between mtDNA and nDNA.
A third explanation would be that a selective sweep of the mtDNA occurred, as already been shown in various animal groups [51][52][53][54][55]. If a mitochondrial selective sweep occurred we should obtain significant results for the neutrality tests only for the mtDNA gene, and a lower nucleotide diversity of the mitochondria compared to the nDNA. This was never observed in our dataset ( Table 4), suggesting that a selective sweep cannot explain the contrasting pattern of sequence divergence observed between mtDNA and nDNA.

Taxonomy
The phylogenetic analysis of Cyt b revealed at least 8 welldefined clades within L. aenigmamus, differing by 5% to 14% K2P sequence divergence. These values are similar to those observed between most mammalian congeneric species [56,57]. According  to Baker and Bradley [57] populations that are distinguished by such genetic distance values merit further study using other molecular loci or methods to determine if an unrecognized species diversity exist. In fact they may represent distinct species under the genetic species concept (a genetic species is defined as a group of genetically compatible interbreeding natural populations that is genetically isolated from other such groups). Our nDNA data also indicate some level of genetic difference between clades. Thus, Laonastes could represent a complex of species. This hypothesis should now be subject to corroboration assessment before we can describe and name them.
Laonastes aenigmamus was described from the vicinity of Ban Muang and Ban Doy, Thakhek District (17.562uN, 104.819uE). Twelve paratypes were also collected near the village of Ban Muang, and two were collected from Thakhek market. Incomplete Cyt b sequences are available in Genbank for three of them: one from Ban Muang (specimen BMNH 1998.25, Genbank AM407933) and two from Thakhek market (BMNH 1998.409 and BMNH 1998.410, Genbank DQ139932 and DQ139933 respectively). When these sequences are included in our phylogenetic analyses, specimen from Ban Muang clusters unambiguously with clade B2, while specimens from Thakhek market unambiguously cluster with clade C. Thus, according to the collection locality of the holotype, the name L. aenigmamus should be retained for specimens of clade B2. Additional morphometric and nuclear data (microsatellites) are needed to definitively conclude on the species-specific status of all other clades. Specimens used in this study are housed at the Museum National d'Histoire Naturelle (Paris, France) and available for additional morphometric and/or molecular analyses (see Table S1 for more details). Multivariate morphometrical analyses of the skulls of all specimens, including the specimens sequenced in this study and the type specimens housed at the BMNH, are needed to investigate intra-and interpopulation morphometrical variability. This will enable us to test morphometrical differentiation between clades and to potentially assign all type specimens to the corresponding clade. Additional molecular data (nuclear DNA sequences, microsatellites) will enable us to further investigate historical and ongoing gene flow between populations. All these data should allow us to evaluate the taxonomic significance of the genetic differentiation observed between clades and to frame it in formal zoological nomenclature.

Geographical distribution
The Khammuan Limestone forms a belt of karsts 290 km long and 30-120 km wide that stretches NW-SE across the full width of Central Laos [58,59]. This formation is isolated from other Table 2. Average percentage of pairwise differences (Kimura 2-parameter) within and between clades (A-H) identified in the phylogenetic analysis, and within and between subclades (A1-A5, B1-B2, D1-D2, E1-E3).  Table 3. Average percentage of pairwise differences (Kimura 2-parameter) within and between subclades A1-A5, B1-B2, D1-D2, E1-E3. limestone, either in northern or southern Laos, by wide plains. To the West, the area is lined by the Mekong River and, to the East, by the Truong Son Mountains, which form the border between Laos and Vietnam. Since 2006, we performed numerous field studies in different provinces in Lao PDR and in Northern Thailand. During this period we questioned villagers about the presence or absence of L. aenigmamus. The Khammuan Limestone area was surveyed with particular attention (Fig. 1). Any time we were moving away from this area, the villager's responses to our questionnaire were becoming progressively negative. When asking people about the presence of Laonastes in small limestone formations peripheral to the central karst area, the answer was that Laonastes was present in the past but is now extinct. Thus, our surveys suggest that L. aenigmamus is only present in Khammuan Limestone area, as described on Fig. 1.

Inferred history of the karsts of the Khammuan Province
Our phylogenetic analyses show that Laonastes populations split into 8 main genetic lineages separated by high genetic distances, and with non-overlapping geographic distributions. Moreover within clades several subclades could be identified with either allopatric (A1 to A5, E1 and E2) or sympatric (D1 and D2) distributions. Moreover E3 and B2 both have a very restricted distribution (captured in only one locality) and were collected in sympatry with E1 and B1, respectively. Figure 4 shows that clades distribution matches karst distribution: an inland insular model for Laonastes population structure is observed, as previously proposed by Rivière-Dobigny et al. [16]. A phenomenon of isolation by distance was detected, but genetic distance may be either low or high between spatially close populations. Thus, isolation by distance, only, cannot explain all the between clade differences. Observations in the field showed that once released in an open space, this species is unable to gallop, due to its remarkable anatomical peculiarities (lumbar vertebrae and pelvis fused together; weak muscles and skeleton) [50]. Because of a low likelihood of survival, if attempting to cross extensive open area, a flat plain of a few hundred meters between two limestone blocks could be sufficient to prevent migration of the animals.
If Laonastes is strictly endemic to karsts of the Khammuan area, its phylogeographical history may be closely tied to the history of its habitat. According to geological studies, Khammuan limestone appeared during the Middle Triassic, and was later buried under Mesozoic sandstone. A new phase of karstification occurred at the beginning of the Tertiary [58][59][60][61]. Aquatic erosion and abrasion by river-borne allogenic clasts contributed to the formation of a particular landscape where limestone pinnacles and towerkarst (up to 300 m high) are prominent and delineate the margins of karstic massifs arising from alluvial plains [58]. Our data strongly suggest that Laonastes current distribution may result from the progressive fragmentation of a formerly panmictic population: Laonastes ancestor being widely distributed within the Khammuan limestone around 8-9 Mya (i.e. at the end of the Miocene). No date for karst fragmentation is given in the bibliography but, if we hypothesize that long-term isolation, lineage sorting and subsequent divergence between Laonastes phylogroups reflect the timing of erosion, as resulting from paleoclimatic events, a possible scenario is as following: -isolation of two south-eastern blocks (corresponding to clade G and H) during the early Pliocene; -separation of a northern block (corresponding to clades E and F) at the end of the Pliocene; -rapid fragmentation in the area corresponding to clades A, B, C and D (which relationships remain largely unresolved), mainly during the Pleistocene.
If the patterns of distribution of this species reflect the general history of the landscape, similar patterns should be discovered in other species with similar ecological constraints and a similar geographical range. The hypothesis of successive karst landscape fragmentations in Khammuan (at the end of the Tertiary or beginning of the Quaternary) may be tested through the biogeographical analysis of other endemics. Of particular interest for this type of study would be the gymnure, Hylomys megalotis and the murid rodent, Saxatilomys paulinae; both known to be endemic to the Khammuan Limestone Province [11,12].

Conservation
Limestone is an important raw material for construction industry and karsts face increased risks of destruction from quarrying activities. As a consequence, extinctions of limestone restricted species resulting of economic development have been recorded [8]. According to Dawson et al. [14] from a paleontological and phylogenetic perspective, efforts to conserve Laonastes, the sole survivor of a morphologically distinctive family of rodents with deep evolutionary roots in Asia, should be given the highest priority. The survival of Laonastes is closely linked to the preservation of its habitat. Moreover our results suggesting that distinct karst towers harbor distinct evolutionary significant unit (that may represent distinct species), it is important to protect not only largest limestone areas but also isolated limestone hills.
Laonastes is not the only vertebrate species known to be endemic to the Khammuan Limestone region where, very probably, many invertebrate and plant species remain to be described and the conservation of Laonastes habitat would also allow the conservation of numerous other species.

Conclusion
Our analyzes of L. aenigmamus molecular data reveal the existence of important inter-populational differences within a restricted area (about 200 km650 km). Limestone karsts are wellknown as ''arks'' of biodiversity and often contain high levels of endemism. However, depending on the dispersal capacities of the taxa, the level of endemicity varies: micro-endemism (i.e. site endemism) is much more frequent in invertebrates than in vertebrates, and in cave-adapted fauna than in surface fauna [4]. Exceptional micro-endemism cases were reported for snails: among just eight selected land snail genera, a large number of species (78) were found to be site endemics (species restricted to single isolated karsts) in Peninsular Malaysia. Concerning mammals a recent study on Leopoldamys neilli showed strong phylogeographic structure with six allopatric genetic lineages corresponding to particular regions of Thailand. However the degree of Cytb (K2P distance) genetic divergence between these clades was much lower than those observed for Laonastes (1.5-7.3%) despite higher geographical distances (20-630 km) and strongly isolated kart blocs [62]. Thus, Laonastes distribution represents an exceptional example of micro-endemism for mammals. Furthermore, the high levels of molecular divergence observed between some populations, suggest that L. aenigmamus may represent a complex of species and/or sub-species.
Laonastes is endemic to the Khammuan Limestone National Biodiversity Conservation Area and is ecologically strongly specialized. It may represent a timeline and spatial reference model to which compare the data collected for other species endemic to this region in order to infer the geological and biological history of the limestone formation in the Khammuan Province. This increases the necessity of a strict protection of this rare animal, which cannot be achieved without the protection of its natural habitat. Figure S1 Phylogenetic relationships between GHR (left) and BFIBR (right) sequences recovered by maximum likelihood (ML) analysis (GTR+G and HKY+G substitution model, respectively). Main clades are identified on the right border. Black dots represent supported nodes (ML Bootsrap.75, Bayesian PP.0.95).

(TIFF)
Table S1 List of specimens used in this study with the field number, mt DNA genetic clades identified in the phylogenetic analyses, GenBank number and locality of origin and GPS coordinates. (XLS)