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

Taxon sampling

For molecular phylogenetic analysis, we included 71 Everettia specimens representing 16 of the 17 known species from Sabah. Besides, five Everettia species from Kalimantan and four Everettia species from Sarawak were also included (Table 1; Figs. 1 and 2). The specimens were obtained from the following depositories: BORNEENSIS at Universiti Malaysia Sabah, the Sabah Parks Museum (SP), Jaap Jan Vermeulen’s private collection (JJ), Leiden, Naturalis Biodiversity Center, Leiden (RMNH, ZMA), the Natural History Museum, London (BMNH), Mohammad Effendi Marzuki’s private collection (ME) and Yansen Chen’s private collection (YSC). Additional materials were obtained under the permits: Sarawak Forestry: NPW.907.4.4 (Jld.14)-31), WL14/2017; and Sabah Parks: TS/PTD/5/4 Jld.54 (112). For an outgroup taxon, we included two specimens of Quantula striata Gray, 1834, which belongs to the sister genus of Everettia within the family Dyakiidae.

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 and 2018 (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³ of foot tissue of single individuals (either fresh, frozen, or kept in ethanol) was extracted with DNeasy^™ 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^™ 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 μl volumes, using 5 μl 10 × reaction buffer (PROMEGA^® or QIAGEN^®), 5 μl two mM dNTP, 6 μl 25 mM MgCl₂, 2 μl for each primer (5 pmol), 26.85 μl 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 directions using the BigDye Terminator Cycle Sequencing Kit v. 3.1 (Applied Biosystems Ltd., Waltham, MA, USA), on an ABI 3100 Genetic Analyser (Applied Biosystems Ltd., Waltham, MA, USA), by Macrogen^® or the BigDye^® Terminator v1.1, v3.0 and v3.1 Sequencing Kit on an Applied Biosystems 3730xl DNA Analyser at MyTACG Biosciences Enterprise.

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) and partition models (Chernomor, Von Haeseler & Minh, 2016) based on the both AIC and BIC that built into IQ-Tree v.1.6.7 (Nguyen et al., 2015; Trifinopoulos et al., 2016). We limited the candidate models to the six models that are available in MrBayes analysis, namely, JC, F81, K80, HKY, SYM and GTR. The results of ModelFinder and partition model suggested different partition schemes and substitution models for respective AIC and BIC selection criteria (Additional File 2). 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²).

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 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.

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 (Manthey et al., 2017).

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, 1997; Gathorne-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, 2011, 2017; Sheldon, Lim & Moyle, 2015).