Molecular phylogenetics and evolutionary history of the endemic land snail genus Everettia in northern Borneo

Borneo has gone through dramatic changes in geology and topography from the early Eocene until the early Pliocene and experienced climatic cycling during the Pleistocene. However, how these changes have shaped the present-day patterns of high diversity and complex distribution are still poorly understood. In this study, we use integrative approaches by estimating phylogenetic relationships, divergence time, and current and past niche suitability for the Bornean endemic land snail genus Everettia to provide additional insight into the evolutionary history of this genus in northern Borneo in the light of the geological vicariance events and climatic fluctuations in the Pleistocene. Our results show that northern Borneo Everettia species belong to two deeply divergent lineages: one contains the species that inhabit high elevation at the central mountain range, while the other contains lowland species. Species diversification in these lineages has taken place before the Pliocene. Climate changes during the Pleistocene did not play a significant role in species diversification but could have shaped contemporary species distribution patterns. Our results also show that the species-rich highland habitats have acted as interglacial refugia for highland species. This study of a relatively sedentary invertebrate supports and enhances the growing understanding of the evolutionary history of Borneo. Species diversification in Everettia is caused by geological vicariance events between the early Miocene and the Pliocene, and the distribution patterns were subsequently determined by climatic fluctuations in the Pleistocene.

history of geological and climatic stability interspersed with periods of upheaval. During the Palaeogene, east and north Borneo was submerged while the rest of Borneo was connected with other parts of Sundaland. Between Eocene and Pliocene, regional tectonic activities have caused the emergence of land and mountain building in Borneo (Hall, 2013), notably: the formation of highlands in central Borneo, the uplifting of Meratus Mountains in southern Borneo, and uplifting of Mount Kinabalu in northern Borneo. The erosion resulting from these mountain-building events have created the land in the northern, eastern and southern parts of Borneo by filling large basins with sediment.
Borneo has been latitudinally stable, and a large part of it has been covered by tropical forest throughout this period (Lumadyo et al., 1993). Widespread evergreen rainforests would have covered much of Sundaland during the early and middle Miocene (De Bruyn et al., 2014). In the Pleistocene, rainforest persisting in some areas of the island were relatively little affected by climatic fluctuations as compared to other parts of Sundaland (Cannon, Morley & Bush, 2009;Wurster et al., 2010;Morley, 2012).
Most of the previous studies use widespread organisms as model taxa to understand how historical processes shaped the genetic and diversity patterns. However, the genetic and diversity patterns of a well-dispersing taxon may be easily diluted and thus impede the interpretations of the events that shape the patterns (Beck & Rüdlinger, 2014;Manthey et al., 2017). Hence taxa that are relatively sedentary and narrowly distributed, and endemic to Borneo are potentially more suitable model organisms. Slow-moving land snails have proven to be excellent model species to understand evolutionary histories at different scales (Davison, 2002;Hugall et al., 2002), which is why we here employ an endemic land snail genus in the context of the historical biogeography of Borneo.
The taxonomy and distribution of the Bornean endemic land snail genus Everettia Godwin-Austen, 1891 has been well documented in northern Borneo (Liew, Schilthuizen & Vermeulen, 2009). This genus is one of the most speciose macro land snails endemic to Borneo. It occupies different habitats from lowland tropical rainforest to highland montane forest, is found in intact forest, degraded forest, swampy forest, coastal forest and islands. A large number of Everettia species are endemic to the highlands of Borneo, and many lowland species show disjunct distribution (Liew, Schilthuizen & Vermeulen, 2009).
In this study, we aimed to provide a temporal framework for the diversification of the different lineages through the reconstruction of a time-calibrated multilocus species tree using relaxed clock models with species distribution modelling during the last glacial maximum (LGM). We examine whether species diversification in Borneo and highland diversity on Mount Kinabalu are either due to Pleistocene climatic fluctuation or earlier Tertiary palaeogeographic events. First, we estimate the phylogeny of Everettia species in Borneo, including species from Sarawak and Kalimantan, by using mitochondrial and nuclear DNA, to elucidate the evolutionary history of Everettia in northern Borneo in the light of the key vicariance events. Second, we construct species distribution models for Everettia species in Sabah, where extensive occurrence data are available, to examine the changes of species distributions during the last glacial period and identify possible refugia during the LGM.
For species distribution modelling, we obtained distribution records of Everettia species from the BORNEENSIS Molluscan collection that consists of 860 collection lots of Everettia species from Sabah that were collected between the years 2000 and2018 (Figs. 3-6). After excluding collection lots for which the exact location and species identity could not be determined, the final distribution data consists of 718 collection lots, which comprise 2,024 specimens of 17 Everettia species from Sabah (Additional File 1). The sampling bias in the distribution data from BORNEENSIS collection is negligible as the entire surface of Sabah has been covered adequately in terms of the geographical space, with some areas having been sampled more densely due to the heterogeneity of the habitat such as mountain ranges and islands (Figs. 3-6).

Molecular methods
Genomic DNA from approximately 2-3 mm 3 of foot tissue of single individuals (either fresh, frozen, or kept in ethanol) was extracted with DNeasy TM nucleic acid extraction kits (QIAGEN Ò , Hilden, Germany) and subsequently stored at −20 C. Then, PCR was performed using a PTC-200 thermocycler (MJ Research, Inc., St. Bruno, QC, Canada) or T100 TM Thermal Cycler (BIO-RAD, Hercules, CA, USA) to amplify the mitochondrial DNA regions 16S with the primer pair 16Sbr-L and 16Sbr-H (Palumbi et al., 1991) and COI with primers LCO1490 and HCO2198 (Folmer et al., 1994). Also, the nuclear rDNA region ITS-1 was amplified with the primer pair 5.8c 'silkworm' and 18d 'fruitfly' (Hillis & Dixon, 1991) and 28S with primers 28S1128 and 28S2119R (De Weerd, 2008). PCR reactions were performed in 50 ml volumes, using 5 ml 10 × reaction buffer (PROMEGA Ò or QIAGEN Ò ), 5 ml two mM dNTP, 6 ml 25 mM MgCl 2 , 2 ml for each primer (5 pmol), 26.85 ml de-ionized autoclaved water and 1 unit of Taq polymerase (PROMEGA Ò or QIAGEN Ò ). Later, the following cycling profile was used: 2 min at 95 C, followed by 35 cycles of 1 min at 95 C, 1 min at 55 C for 16S, COI and 28S (60 C for ITS-1) and 2 min at 72 C, and a final extension period of 10 min at 72 C. Next, PCR-amplified DNA fragments were purified with the High Pure PCR Product Purification Kit (Roche Ò or ExoSAP-IT Ò ), according to the manufacturer's protocol. Finally, DNA sequencing was performed directly on purified PCR products in both

Phylogenetic analysis
A total of 96 genetic sequences of the previous study (Liew, Schilthuizen & Vermeulen, 2009) and 160 new genetic sequences from the present study were aligned using the ClustalW multiple alignment algorithm in the BioEdit Sequence Alignment Editor, version 7.0 (Hall, 1999) and manually adjusted with the same programme. Before the phylogenetic analyses, the data matrix was partitioned by markers and codons of COI, namely, first, second and third codon positions of COI, 16S rDNA, ITS-1 and 28S rDNA. Then, each of the partitions was tested for molecular evolution via ModelFinder (Kalyaanamoorthy et al., 2017)  We explored the phylogenies estimated based on different substitution models selected for AIC and BIC but the resulted phylogenies are generally congruent (Additional File 3). Hence, we used the best-fit substitution models and partition scheme of BIC selection: partition (1) 16S+ITS: GTR+F+G4, partition (2) COI1+COI2+28S: SYM+I+G4 and partition (3) COI3: GTR+F+G4. The sequences were analysed using Bayesian analysis (BA) with MrBayes 3.1 (Huelsenbeck & Ronquist, 2001) at the CIPRES Science Gateway portal (Miller, Pfeiffer & Schwartz, 2010) and a maximum likelihood (ML) method implemented in IQ-Tree v.1.6.7 (Nguyen et al., 2015). For BA, the data matrix was analysed with 10 million generations and sampled every 1,000th generation. Then, we discarded the first 25% of the samples. BA was repeated three times for data matrix, and a consensus tree with a cut-off value of 50% was calculated for the resultant trees. For ML analysis, we estimated the phylogeny by using 1,000 ultrafast bootstrap replicates (Minh, Nguyen & von Haeseler, 2013).

Estimation of divergence time
BEAST 2 (ver. 2.6.1) (Drummond & Rambaut, 2007) was used to estimate the timescale for Everettia species divergences based on selected samples for each species. We presume that the split between two Everettia species: E. sp. 1 and E. sp. 2, that occur at the two sides of the Meratus range in South Kalimantan based on a geological event -the uplift of the Meratus Range during late Miocene (10 Ma) (Hall, 2013). Hence, the hypothesis on the timing of speciation of the phylogeny is based on this calibration point which the divergence of the species has resulted from the uplifting of the mountain ranges in Borneo. The tools provided in BEAST 2 were used to estimate node ages to the most common recent ancestor of the split and substitution rates. We carried out four independent runs of 50,000,000 generations each, sampled every 10,000 generations, using calibrated Birth-Death model with best-fit GTR models, a relaxed lognormal molecular clock was employed, and default options for all other priors and operator settings. The Birth-Death model is chosen as we believe that the evolution of Everettia species a continuous-time process with a probability that a lineage will go extinct. We also explored the time divergence estimates for the combinations two different best-fit substitution models (selected by BIC and AIC criteria) and two calibrated models (Yule model vs Birth-Death model) and the results of these analyses are similar (Additional File 4). The output of each independent run was visualised in Tracer 1.4. Samples and trees from separate runs were pooled after removing the first 10% as burn-in using LogCombiner ver. 2.6.1 and 10% of the trees were discarded as burn-in, and maximum clade credibility trees were calculated each from the remaining 180,004 trees using TreeAnnotator 2.6.0. Divergence dates were computed using BEAST 2 at CIPRESS. The geology-based calibration point (10.0 Ma ± 0.5, 95% CI) was taken as the central trend of a normally distributed prior in BEAUti.

Ecological-niche modelling
To understand how the distribution of Sabah Everettia species has changed during the paleoclimatic fluctuations in the Pleistocene, we predicted ecological niches for all eighteen Sabah Everettia species by using current distribution data under the contemporary (i.e. interglacial) and past (i.e. glacial) climatic conditions. As in other land snail studies (Hugall et al., 2002), we assumed niche conservatism for Everettia.
For the environmental data, we used the bioclimatic dataset version 1.4 (http://www. worldclim.org/current; Fick & Hijmans, 2017). Each of the current bioclimatic layers of resolution of 30 arc-s was clipped to the extent of Borneo. After that, we sampled bioclimatic variables for 500 random locations in Borneo to evaluate the collinearity among the 19 climatic variables by using pairwise Pearson's r correlation (Additional File 5). After we removed highly correlated variables (r > 0.8), a total of seven climatic variables were used for species distribution modelling, namely, BIO1 Annual Mean Temperature, BIO3 Isothermality, BIO4 Temperature Seasonality, BIO7 Temperature Annual Range, BIO12 Annual Precipitation, BIO15 Precipitation Seasonality and BIO19 Precipitation of Coldest Quarter. Next, the corresponding seven bioclimatic variables of the paleoclimatic dataset for the LGM (model CCSM; http://www.ccsm.ucar.edu/, Kiehl & Gent, 2004) were resampled at resolutions of 30 arc-s (~1 km 2 ).
Then, MaxEnt software (ver. 3.4.1, Phillips, Anderson & Schapire, 2006;Phillips & Dudík, 2008) was used to generate logistic probability maps of species presence with logistic values ranging from 0 (unsuitable) to 1 (optimal habitat). The model was run using the following settings: the maximum number of background points = 10,000; replicates = 10; and replicate run type-Cross validate. All other parameters were kept at default values. Finally, the average of the logistic probability of species occurrence for each grid cell was calculated from the resultant ten replicates.

Phylogenetic analyses
The combined mitochondrial and nuclear DNA matrix comprises 73 specimens and 2,795 characters (16S: 1-501; COI: 502-1059; 28S: 1060-1869; ITS: 1870-2795 (Additional File 6). The best nucleotide substitution models are reported in Additional File 2. As revealed by the Bayesian posterior probability (PP) and maximum likelihood analysis bootstrap (BS) values of the phylogenetic tree in Fig. 7, most of the species are monophyletic, and phylogenetic relationships between species are similar to those found in a previous study (Liew, Schilthuizen & Vermeulen, 2009).
In contrast to the previous study (Liew, Schilthuizen & Vermeulen, 2009), this study shows the phylogenetic relationship of Sabah Everettia species in the broader context of Bornean Everettia species. Everettia species of Sabah do not form a monophyletic group, and belong to four independent lineages, namely: lineages A, C, D and E (Fig. 7). The other lineage B consists of one species from Sarawak near the border with Brunei and two species from South Kalimantan. However, some of the phylogenetic relationships among these lineages are poorly supported by Bayesian analysis (i.e. PP < 0.95) (Fig. 7).
A total of 12 out of 16 Everettia species in Sabah belong to two major lineages. The first lineage (hereafter, lineage A) consists of nine species, seven of which are lowland species that have their lowest elevation distribution below 1,000 m, namely E. subconsul, E. interior, E. paulbasintali and E. jucundior from Sabah (Figs. 5-7); E. algaia from Sarawak; E. sp. 4 and E. sp. 5 from East Kalimantan (Fig. 1). Two of the species of this lineage (E. layanglayang and E. themis) have their lowest elevational distribution below 2000 m (Figs. 4 and 7). With this expanded genetic dataset, E. themis is now paraphyletic to E. subconsul.
The second lineage (hereafter, lineage C) consists of eight species, of which four are Mount Kinabalu endemics with a lowest elevational limit above 2,000 m, namely E. jasilini, E. safriei, E. corrugata corrugata, and E. c. williamsi (Figs. 3 and 7); two are highland species with their lowest elevation above 1,000 m, namely E. monticola and E. dominiki; and a further two are lowland species: E. planispira from Sabah and E. sp. 3 from Central Kalimantan which occur more than 600 km apart from each other (Figs. 1 and 6).
The remaining four Sabah Everettia species, namely E. jucunda, E. klemmantanica, E. lapidini and E. jucundior, do not belong to the lineages A and C. The Sabah Everettia jucunda form a lineage with an Everettia species (sp. 6) from Sarawak. The Sabah and Sarawak species are more than 500 km apart from each other (lineage E, Fig. 7). E. lapidini and E. klemmantanica are not shown as mutually monophyletic species but as a joint monophyletic clade (lineage D, Fig. 7).
The lineage B consists of two Everettia species from South Kalimantan (sp. 1 and sp. 2) and E. baramensis from Sarawak. The Sarawak and South Kalimantan species are more than 700 km apart from each other (lineage B, Fig. 7). Lastly, E. sp. 7 from Sarawak does not form a clade with any other Everettia species.

Divergence time and tempo of speciation
Here, we used only one calibration point based on a single biogeographic event given the limited availability of fossil records, and reasonable estimates of mutation rates across Figure 7 The phylogeny of 25 Everettia species with Quantula striata as outgroup. Bayesian inference 50% majority-rule consensus trees based on the concatenated dataset consisting of parts of 28S, ITS-1, COI and 16S. Bayesian posterior probabilities and bootstrap support after 1,000 maximum likelihood replicates are shown above and below the branches of the nodes. The font and colour of the taxa name on the tree indicate the distribution of the species. The colour panels next to the taxa names indicated the lowest elevation distributional (Left) and highest elevation distributional (Right) of the species. The number after the taxa name specimen number of different genes for the land snail taxa in this region. Currently, the only known land snail fossils in Southeast Asia are from species of Family Cyclophoridae that cannot be used for calibration in this study (Raheem et al., 2018;Xing et al., 2019). The topography of the chronogram is generally congruent with the phylogenetic analysis, of which most of the deeper nodes are poorly supported (PP < 0.95) ( Fig. 7; Fig. 8). Our results show the divergences among Everettia in various areas of Borneo are tally to the area's major mountain uplifting events. These divergence time estimates are based on the hypothesis that mountain uplifting events caused the divergence of the two Everettia species at the two sides of the Meratus range could be falsified in the future if there are more accurate vicariance geological events or reliable fossil record available to improve the calibration of the phylogeny. Diversification of Everettia species in Borneo began in the Late Oligocene (25.8 Ma). These species diversified into five major lineages between the early Miocene (23-17 Ma). The lowland lineage (lineage A) diversified rapidly into seven species between 7 and 19 Ma (Fig. 8). The highland lineage (lineage C) diversified rapidly into montane species and the Mount Kinabalu endemics lineage between 4 and 15 Ma (Fig. 8). Deep divergence of the South Kalimantan and Sarawak species is seen in lineage B (17 Ma).

Comparison of the ecological-niche model of contemporary and past distribution
Distributional range shifts of Everettia species during the LGM are predicted by the habitat suitability maps in Figs. 3-6. As shown in the phylogenetic analysis, E. themis is now considered as E. subconsul for species distribution modelling (Fig. 7). The area under the curve (AUC) values for the 15 species models are higher than 0.85, except E. klemmantanica.
Most of the Sabah Everettia species have their suitable habitat in Sabah, particularly the endemic species on Mount Kinabalu and central mountain ranges in Sabah. The analysis suggests that suitable habitats for E. jucundior, E. planispira and E. interior are not limited to Sabah, but are extended to large areas in the eastern and southern part of Borneo. Besides, small areas of suitable habitats for E. paulbasintali, and E. occidentalis are located in the eastern part of Borneo.
The palaeoclimatic models predict contraction and expansion of suitable habitats during the LGM for different Everettia species. All four Mount Kinabalu endemic species (Everettia corrugata corrugata, E. c. williamsi, E. jasilini and E. safriei) have experienced range expansion during the LGM at the central mountain range of Sabah. Highland species E. dominiki, E. monticola, E. layanglayang and E. lapidini experienced range expansion as Mount Kinabalu endemics and also in the mountain ranges in the western Borneo.
Everettia planispira, the lowland relative of E. dominiki and E. monticola-experienced significant range expansion in eastern and southern Borneo. Phylogenetic analysis suggests that E. planispira is the sister taxon for E. sp. 3, which is found in southern Borneo. A few of the lowland species, viz. E. paulbasintali, E. occidentalis, E. jucunda, and E. jucundior, experience range contraction and probably remain with very limited suitable habitats. E. subconsul was predicted to have experienced a shrinking of suitable habitat during the LGM into areas near the tip of northern Borneo, including offshore islands and lowland around Mount Kinabalu. Phylogenetic analysis also showed that the populations of E. subconsul on northern offshore islands and tips of northern Borneo are the oldest for the species. The other lowland species E. interior experienced a little reduction of suitable habitats and its contemporary distribution range is similar to that during the LGM. In particular, the contemporary distributional range of E. interior could potentially extend to eastern Borneo.

DISCUSSION
A high species diversity and high degree of endemism in northern Borneo are well known for many plant and animal taxa, particularly for the central mountain ranges, that is the Crocker Range, Mount Kinabalu and the Trusmadi Range (Liew, Schilthuizen & Vermeulen, 2009;Beaman, 2005). Land snail studies in other regions suggest that vicariance events that persist long enough play crucial roles in driving radiation (Douris et al., 1998;Parmakelis et al., 2005;Fiorentino et al., 2010;Pfenninger et al., 2010;Rowson, Tattersfield & Symondson, 2011), with other factors such as dispersal events and niche differentiation causing further modification (Douris et al., 1998;Schilthuizen et al., 2004;Hausdorf & Hennig, 2004Holland & Cowie, 2009;Ketmaier et al., 2010;Kokshoorn et al., 2010). Previously, the phylogeny of Everettia species was estimated without other congener species from outside of Sabah (Liew, Schilthuizen & Vermeulen, 2009). Although species sampling outside of Sabah is still far from complete, these additional species from part of Borneo provides a more accurate phylogeny to illustrate the evolution of Sabah Everettia species that are more or less completely sampled reveal several novel insights.

Divergence of species in the highland lineage
First, most of the Sabah species belong to two deeply diverged lineages. One lineage mainly consists of highland species, particularly all endemics of Mount Kinabalu, while the other lineage includes lowland species. The divergence of these two lineages took place during the early Miocene, which coincided with the uplift of mountain ranges and an extended land area from the southwest to the northeast of the centre of Borneo (Fig. 8). Hence, the divergence was not caused by the more recent uplift of Mount Kinabalu as postulated by studies on other organisms (O'Connell et al., 2018).
The diversification of the four Kinabalu endemics (E. jasilini, E. safriei, E. corrugata and E. c. williamsi) within the highland lineage happened after the middle Pliocene (after 3.8 Ma), and could have been caused by the uplifting of Mount Kinabalu (Figs. 3 and 4). The rapid uplift of Mount Kinabalu at the rate of 500 m per million years (Cottam et al., 2010) could have caused allopatric speciation when the habitat at higher elevation arose, and populations were isolated (Merckx et al., 2015). However, the remaining three species (E. monticola, E. dominiki and E. layanglayang) that reach to an elevation of 3,000 m on Mount Kinabalu and are sympatric with the four Kinabalu endemics more likely diverged by geographical isolation on other mountain summits and subsequently became secondarily sympatric (judged by their deep divergence, before the emergence of Mount Kinabalu).
The palaeo-distributions during the LGM of these seven species provide some insights that these species had more widespread distribution ranges in the central mountain ranges of Borneo that are adjacent to Sabah, based on the suitable habitat analysis (Figs. 3 and 4). This suitable habitat may have facilitated dispersal of these once geographically isolated highland species between Central and northern Borneo montane areas when the cooler temperature during the LGM caused the montane forest to descend and spread, which would have increased connectivity among mountains .
However, habitat at lower elevations became hostile to these highland lineage species when the climate warmed up during interglacials. These species probably reacted by moving to suitable habitat at higher altitudes or went extinct altogether. Thus, we believe that Mount Kinabalu has served as a refugium during interglacial periods for highland Everettia species. These highland species could have been trapped there during several glaciation cycles, although we cannot say at which Quaternary glaciation stages this happened. Furthermore, we have shown that land snails on other northern Bornean mountains also show shorter ranges at higher elevations compared to the lowland and lower montane areas (Liew, Schilthuizen & Lakim, 2010), indicating that these species have been pushed upwards until the end of their optimum habitat. This finding supports the studies of other taxa that proposed the mountain ranges in Sabah play a role in the maintainance of ancient lineages (Sheldon, 2017).
The discrepancy of the two divergent processes for the sympatric species on top of Mount Kinabalu provides additional insight that challenge the conventional view that Mount Kinabalu acted as a 'speciation pump' and that lower elevation ancestors gave rise to high-elevation endemics (Lee & Lowry, 1980;De Laubenfels, 1988;Holloway, 1996;Chan & Barkman, 1997;Barkman & Simpson, 2001;Tanaka et al., 2001).

Divergence of species in the lowland lineage
Sabah became fully emergent only at the end of the Miocene or Early Pliocene. Two of the most widespread lowland species in Sabah-E. subconsul and E. paulbasitali from lineage A, rapidly colonised newly emerged habitat. Although we did not perform analysis on different populations of other lowland species, we think it is very likely the other widespread lowland species, for example E. jucundior and E. planispira, dispersed to the newly formed land at the same time. In addition to the role of Mount Kinabalu as an interglacial refugium for highland lineage species, SDM analysis shows that Mount Kinabalu also acted as a glacial refugium for lowland lineage species, for example E. subconsul. northern Borneo has been mentioned as a probable glacial refugium during climate changes in the Pleistocene (Brandon-Jones, 1996, 1997Gathorne-Hardy et al., 2002), but the exact locations of suitable refugia have remained unknown, with some hypothesising that Mount Kinabalu and the Crocker Range could have played such a role (Cockburn, 1978;Smith, 1980;Quek et al., 2007;Jalil et al., 2008). Our study identifies two probable glacial refugia for E. subconsul, on the east and west slopes of Mount Kinabalu. These two glacial refugia, together with unsuitable habitat and mountain ranges as geographical barriers in the centre of northern Borneo, could explain how the east and west coast populations of E. subconsul have maintained their deeply diverged origin since the late Miocene (Figs. 6 and 8).
In contrast to the distribution patterns in the highland lineage, most of the species in lowland lineage occur allopatrically, with the exception of E. subconsul, E. paulbasintali, E. planispira and E. jucundior which are sympatric on the east coast of Sabah. At first glance, the allopatric distributions of the lowland Everettia species appear to be due to geographical isolation caused by mountain ranges, as has been suggested in studies on other taxa (Bänfer et al., 2006). Besides, distribution patterns of lowland species are similar to physiography, vegetation and biozoographical subregions of northern Borneo (Collenette, 1963;Mackinnon & Mackinnon, 1986;Mackinnon et al., 1996;Wong, 1998).
Based on the palaeo-distribution analysis, the lowland species mostly expanded post-glacially, whereas the ranges of the highland species are currently contracting by moving to higher elevations. These different responses by highland and lowland land snails to climate fluctuations are also known from other tropical regions (Wronski & Hausdorf, 2008). Hence, Mount Kinabalu acts as interglacial and glacial refugium for remnant populations, which results in a species diversity hotspot. In Everettia, a total of thirteen out of seventeen northern Borneo species occur on Mount Kinabalu, and six of those species are endemic. The high richness of ancient species agrees with the fact that northern Borneo has had a stable ever-wet climate with most of the forest persisting over the glaciations (Bird, Taylor & Hunt, 2005;Wurster et al., 2010). northern Bornean populations or taxa are known to have been isolated from other parts of Borneo, especially western Borneo, in rainforest refugia during the Pleistocene (Moyle et al., 2005(Moyle et al., , 2011Sheldon, Lim & Moyle, 2015).

CONCLUSIONS
Our data enhance the understanding of the evolutionary history of northern Borneo. The northern Borneo Everettia species belong to two deeply diverged lineages. The ecological differentiation and divergence of these two lineages were caused by the uplift of mountain ranges in central Borneo during the Miocene. The continuing eastward and northward extension of Borneo land area together with the formation of central mountain ranges in these newly emerged parts of Borneo have probably driven the species diversification of Everettia in both lineages throughout the Miocene. The species distributional ranges have changed during fluctuating climatic conditions in the Pleistocene. The highland species tended to expand their distribution ranges and lowland species distributional ranges retracted in response to glacial periods, and vice versa during interglacials. We also show that the central mountain ranges of northern Borneo, especially Mount Kinabalu, have acted as refugia in both interglacial and glacial periods. Thus, the contemporary species richness and endemism are caused by geological vicariance events while the contemporary species diversity and distribution patterns are shaped by the Pleistocene climatic fluctuations. We also provide a scenario for how these mountain ranges may have served as refugia for lowland and highland species during both warm and cooler periods. In fact, less than 1% of the total land surface of Borneo is above 2,000 m, and more than three-quarters of this is in northern Borneo. Hence, highland habitats are importance as future refugia for species impacted by rapid climate change in the near future.