Total genomic DNA was isolated from liver samples using conventional methods. The protocol was adapted from Yari et al. (2017). A 1494 bp fragment of mitochondrial gene (mtDNA) NADH dehydrogenase subunit 2 (ND2), including the flanking tRNAs, and 1077 bp fragment of nuclear gene (nuDNA) recombination activating gene 1 (RAG1) were amplified through PCR. The primers used for amplification and sequencing are listed in Table 1. The PCR solution was prepared by adding 2 μl of genomic DNA, 2 μl of forward primer, 2 μl of reverse primer, 25 μl of NEXpro e PCR 2X Master Mix, and 19 μl of ddH₂O into a 0.5 ml microtube that labelled with an identification number. The PCR reactions for ND2 were performed on a Biometra TAdvanced thermal cycler under the following conditions: initial denaturation at 95˚C for 2 min, followed by a second denaturation at 95˚C for 35 s, annealing at 52˚C for 35 s and a cycle extension at 72˚C for 35 s, for 31 cycles, and a final extension at 72˚C for 5 min. The PCR reactions for RAG1 were performed under following conditions: initial denaturation at 95˚C for 2 min, followed by a second denaturation at 95˚C for 35 s, annealing at 50˚C for 35 s and a cycle extension at 72˚C for 95 s, with four seconds added to the extension per cycle, for 34 cycles, and a final extension at 72˚C for 5 min. The PCR products were visualised on 1% agarose gel electrophoresis and sent to Apical Scientific (Malaysia) for purification and sequencing. The sequences were analysed and edited in MEGA v11 (Tamura et al. 2021).

Besides the specimens collected in the field, additional sequences used in the analyses were obtained from GenBank (Table 2). Hemiphyllodactylus sequences were either selected directly from the BLAST results or manually searched by species name or their published GenBank accession number. Outgroups were selected by following Grismer et al. (2015a). The sequences were then aligned by MUSCLE (Edgar 2004) in MEGA v11 (Tamura et al. 2021). Collection abbreviations as follows: ABTC = Australian Biological Tissue Collection, Australia; ACD = Arvin C. Diesmos field collection; AMB = Aaron Bauer field series; AMS = Australian Museum, Sydney, Australia; BNHS = Bombay Natural History Society, Mumbai, India; CAS = California Academy of Sciences, USA; CBC = Centre for Biodiversity Conservation, Royal University of Phnom Penh, Cambodia; CCA = Louisiana State University Museum of Natural Science, USA; CES G = Karanth lab field series and CES L = Kartik Shanker lab field series at the Centre for Ecological Sciences, Bangalore, India; CSUFT = Central South University of Forestry and Technology, Changsha, China; DWB = Donald W. Buden; FK = Fred Kraus field series; FMNH = Field Museum of Natural History, USA; ID = Indraneil Das field series; IEBR = Institute of Ecology and Biological Resources, Hanoi, Vietnam; ITBCZ = Institute of Tropical Biology Collection of Zoology, Vietnam; JAM = Jimmy A. McGuire field series; JB = Jon Boone; JEM = John E. Measey; KIZ = Kunming Institute of Zoology, China; LSUHC = La Sierra University Herpetological Collection, USA; MCZ = Museum of Comparative Zoology, Harvard University, USA; MVZ = Museum of Vertebrate Zoology, Berkeley, USA; NCBS-BH = National Centre for Biological Sciences, Bangalore, India; NCSM = North Carolina Museum of Natural Sciences, Raleigh, USA; NJNUh = Nanjing Normal University, Jiangsu, China; NUM = Nagoya University Museum, Japan; NUOL = National University of Laos, Laos; RAH = Rod Hitchmough; ROM = the Royal Ontario Museum, Canada; SAM = South Australian Museum, Australia; SYS = Museum of Biology, Sun Tat-sen University, Guangzhou, China; TG = Tony Gamble; TNHC = Texas Natural History of Collection, USA; USMHC = Universiti Sains Malaysia Herpetological Collection, Malaysia; USNM = United States National Museum, USA; USNM-FS = United States National Museum, Field Series, USA; VNUF = Vietnam National University of Forestry, Vietnam; USNM = Zoological Museum, Amsterdam, Netherlands; ZMKU = Zoological Museum, Kasetsart University, Thailand; ZMMU = Zoo­logical Museum of Moscow University, Russia; ZRC = Zoological Reference Collection, Raffles Museum of Biodiversity Research, Singapore; ZSM = Zoologische Staatssammlung, München, Germany.

Phylogenetic trees of total 2571 bp of concatenated mtDNA and nuDNA datasets were constructed using Maximum Likelihood (ML) and Bayesian Inference (BI). In addition, phylogenetic trees using only the nuDNA dataset was constructed in order to test Cobos et al.’s (2016) hypothesis. The best-fit models of evolution for each codon position (Table 3) were calculated using Akaike Information Criterion (AIC) scores recovered in Partition Finder v2 (Lanfear et al. 2017). Maximum Likelihood (ML) analysis was performed using IQ-TREE v2 (Minh et al. 2020) via command line. The nodal supports were estimated by employing 1000 bootstrap pseudoreplicates via the ultrafast bootstrap (UFB) approximation algorithm (Hoang et al. 2018). Bayesian Inference (BI) analysis was performed in MrBayes v3.2 on XSEDE (Ronquist et al. 2012) using CIPRES (Cyberinfrastructure for Phylogenetic Research; Miller et al. 2010). Two simultaneous runs were performed with four chains per run, three hot and one cold using default priors. The analysis was run for 100 million generations and sampled every 10,000 generations from the Markov Chain Monte Carlo (MCMC). Their convergence was verified in Tracer v1.7 (Rambaut et al. 2018) by checking that the effective sample sizes (ESS) were all above 200. The first 10% of the trees from each run were discarded as burn-in. A consensus tree was computed from the two parallel runs using the sumt function in MrBayes v3.2 (Ronquist et al. 2012). Nodes with UFB values of ≥95 and BI posterior probabilities (BPP) of ≥0.95 were considered strongly supported, while UFB and BPP from 90–94 and 0.90–0.94 were considered well-supported (Huelsenbeck and Ronquist 2001; Wilcox et al. 2002; Minh et al. 2013). Uncorrected pairwise sequence divergences were calculated in MEGA v11 (Tamura et al. 2021) using the pairwise deletion option.

A fossil-calibrated chronogram of the concatenated mtDNA and nuDNA datasets was constructed using BEAST v2.7 using an input file generated in BEAUTi v2.7 (Bouckaert et al. 2019). The models of molecular evolution selected are listed in Table 3. The fossilised birth-death model was applied as it yields more accurate node age estimates and does not become computationally expensive even when many calibrations are applied (Heath et al. 2014). A lognormal relaxed clock and the fossil calibrations were applied following Grismer et al. (2015a). The divergence between Sphaerodactylus roosevelti Grant, 1931 and S. torrei Barbour, 1914 was calibrated (exponential, mean = 3, offset = 15) from an amber-preserved fossil Sphaerodactylus from Hispanola dated to 15–20 Mya (Iturralde-Vinent and MacPhee 1996). The divergence between Oedura and Woodworthia was calibrated (exponential, mean = 17, offset = 16) from a fossil from New Zealand ‘Hoplodactylus’ which dates to 16–19 Mya (Lee et al. 2009). The divergence between Pygopus and Lialis was calibrated (exponential, mean = 10, offset = 20) from the fossil Pygopus which dated to 16–19 Mya (Hutchinson 1998). An additional calibration was applied for Gekkota (exponential, mean = 5, offset = 99) following Agarwal et al. (2019), based on amber fossils from Myanmar (Daza et al. 2016). The analysis was run for 100 million generations and sampled every 10,000 generations from the MCMC. The analysis was performed using BEAST v2.7.7 on XSEDE (Ronquist et al. 2012) in CIPRES (Miller et al. 2010). Convergence was verified by noting the ESSs were all above 200 using Tracer v1.7 (Rambaut et al. 2018). The first 10% of the trees were discarded as burn-in. A maximum clade credibility tree (MCC) using mean heights was constructed using TreeAnotator v2.7 (Bouckaert et al. 2019).

The time-calibrated BEAST tree was converted to newick format and pruned to include only the harterti group using the drop.tip function in the ape package (Paradis et al. 2024) in R studio v4.3.3, with one individual per species or population representing each tip. A species occurrence geography file was constructed utilising four discrete areas: the Banjaran Timur located in the eastern one-half of Peninsular Malaysia; the Banjaran Titiwangsa, the central mountain range of the country; the Banjaran Bintang, the northwesternmost range in the country; and the isolated mountain tops Gunung Raya on Langkawi Island, Penang Hill on Penang Island and Gunung Jerai located on the northwestern coast of Peninsular Malaysia between Langkawi and Penang islands. Although Penang Hill and Gunung Jerai are separated by the Straits of Malacca, they were considered a single geographic unit because they are part of the same submerged low-elevation mountain range which is separate from the adjacent Banjaran Bintang to the east. These areas are visible in the relief map in Figure 1. Each species or population was allowed to occupy a maximum of two areas. The ancestral range at each node was estimated using the BioGeoBEARS package (Matzke 2018) in R studio v4.3.3. There are six biogeographical inference models available in BioGeoBEARS, which compose the three main models Dispersal-Extinction Cladogenesis (DEC) model of LAGRANGE (Ree and Smith 2008), Dispersal Vicariance Analysis DIVA (DIVALIKE) (Ronquist 1997) and Bayesian-based BayArea (BAYAREALIKE) (Landis et al. 2013), and each with the option to incorporate founder-effect speciation (+J). The model with lowest AIC and AICc score, and with highest likelihood (LnL) score was selected as the best-fit model.

Morphometric and meristic data were taken, and additional diagnostic characters were recorded following Zug (2010) and Grismer et al. (2020). Morphometric data were taken using Mitutoyo dial callipers to the nearest 0.1 mm on the left side of the body where appropriate. These include snout-vent length (SVL, measured from the tip of snout to the vent), trunk or axillia-groin length (AG, measured from the posterior margin of the forelimb at its insertion point of the body to the anterior margin of the hindlimb at its insertion point on the body), head length (HL, measured from the posterior margin of the retroarticular process of the lower jaw to the tip of the snout), head width (HW, measured at the angle of the jaws), eye diameter (ED, the greatest horizontal diameter of the eyeball), snout-eye length (SN, measured from the anteriormost margin of the eyeball to the tip of snout), nares-eye length (NE, measured from posterior edge of the external nares to the anterior margin of the eyeball), and snout width (SW, measured between the external nares).

Meristic characters taken include circumnasal scales (CN, number of scales contacting the nares), postrostrals or intersupranasals (IS, number of scales between the supranasals), chin scales (CS, number of scales contacting the internal edge of the infralabials and mental from the juncture of the second and third infralabials on both sides), supralabial scales (SL, counted from the largest scale immediately posterior to the dorsal inflection of the posterior portion of the upper jaw to the rostral scale), infralabial scales (IL, counted from the largest scale immediately posterior to the dorsal inflection of the posterior portion of the upper jaw to the mental scale), dorsal scales (DS, number of longitudinal dorsal scales at midbody contained within one eyeball diameter), ventral scales (VS, number of longitudinal ventral scales at midbody contained within one eyeball diameter), lamellae beneath the first finger (FL1, number of subdigital lamellae wider than long on the first finger), lamellae beneath the first toe (TL1, number of subdigital lamellae wider than long on the first toe), lamellar formulae (number of U-shaped, subdigital lamellae (split and single) on the digital pads of digits II–V of the hand and feet), femoroprecloacal pores (FP, number of the contiguous rows of femoral and precloacal pore-bearing scales), and the number of cloacal spurs (CloacS) on the side of the hemipenial swelling or immediately posterior to the vent in females.

Colour pattern characters evaluated include the presence or absence of dark pigmentation in the gonadal tracts and caecum, configuration of the dark markings on the dorsum (dorsal pattern), presence or absence of a dark pre- and / or postorbital stripe extending to at least the neck, presence or absence of a linear series of white postorbital and dorsolateral spots on trunk, presence or absence of dark ventrolateral stripe on trunk and the presence or absence of light-coloured, anteriorly projecting arms of the light-coloured postsacral marking.

The morphological data for other species of the harterti group were obtained from Cobos et al. (2016) and Grismer et al. (2013, 2015a), except for H. titiwangsaensis. Hemiphyllodactylus titiwangsaensis is a species complex comprised of three different populations—Cameron Highlands, Fraser’s Hill, and Genting Highlands. As the data from previous studies might be a mixture of all three populations taken by different researchers, they were excluded and the morphological characters for each population of the H. titiwangsaensis complex were recorded and analysed separately. As the specimens were deposited in different collections (LSUHC and USMHC), the measurements for specimens in each collection were taken by different observers. The morphological characters of the newly discovered populations were compared with other species within the harterti group to search for discrete and overlapping diagnostic characters. Additional multivariate analyses were conducted for the populations from Gunung Raya, Langkawi Island; Gunung Jerai; and their closest relative, H. cicak.

Prior to the subsequent analyses, morphometric characters were size-corrected in order to remove any potential effects of allometry, using the package GroupStruct (Chan and Grismer 2022). For characters with left and right values, only the left values were used for the analyses as the right values for H. cicak are unavailable in Cobos et al. (2016). A Levene’s test was conducted to assess homogeneity of variances using the package car (Fox and Weisberg 2019), and a Shapiro-Wilk test was conducted to assess normality. An analysis of variance (ANOVA) was conducted to search for characters with statistically significant mean differences (p value ≤ 0.05) among the populations. Normally distributed (p value > 0.05) characters with equal variances (p value > 0.05) were subjected to a classic one-way ANOVA, while those with unequal variances (p value ≤ 0.05) were subjected to a Welch’s ANOVA, both using the R studio package userfriendlyscience (Peters 2017). Characters not normally distributed (p value ≤ 0.05) were subjected to a Kruskal-Wallis H test, regardless of the equality of variance. Following that, post hoc tests were conducted on characters with statistically different means (p value ≤ 0.05) to identify which population pairs differed significantly from each other for which character. The post hoc test following ANOVA is Tukey’s Honestly Significant Difference (HSD) test, and that following Welch’s ANOVA is Games-Howell test, both of which were conducted using the package userfriendlyscience (Peters 2017). The post hoc test following Kruskal-Wallis test is Dunn’s test, which was conducted using the package dune.test (Dinno 2024). Boxplots and violin plots were then generated for meristic and size-corrected morphometric characters, respectively, in order to visualise the range, mean, and degree of differences between these pairs. All analyses were conducted in R studio v4.3.3.

Multiple factor analysis (MFA) (Grismer 2025) was conducted to visualise the morphospatial positions of each population relative to one another and the clustering of the sampled individuals, using the package FactorMineR (Husson et al. 2017). Femoroprecloacal pore counts were excluded from the analysis as they are only present in males. Ten meristic (CS, CN, IS, SL, IL, DS, VS, FL1, TL1 and CloacS) and six morphometric characters (SVL, AG, HL, SN, HW, and ED) were used for the MFA. Categorical characters were not included as they are either unavailable in Cobos et al. (2016) for H. cicak, or score the same across all individuals of the Langkawi and Gunung Jerai populations and H. cicak. See File S1: MFA_dataframe for the data frame used to perform MFA. The results of MFA were then visualised using the package Factoextra (Kassambara and Mundt 2020). A non-parametric permutational multivariate ANOVA (PERMANOVA) using the package vegan (Oksanen et al. 2024) was conducted to determine if the centroid locations of each population were statistically different from one another. All analyses were conducted in R studio v4.3.3.

The results of the concatenated mtDNA and nuDNA datasets show that the two newly discovered populations from Langkawi Island and Gunung Jerai are embedded within the harterti group (Fig. 2). Besides that, an additional population (H. cf. titiwangsaensis) that highly resembles H. titiwangsaensis morphologically, but genetically distinct was recovered. The ML and BI analyses produced phylograms with generally similar topologies with strongly supported nodes, except for H. harterti + (H. cicak + (Langkawi population + Gunung Jerai population)) (BPP/UFB: 0.72/89), and the position of H. cf. titiwangsaensis) that highly resembles H. titiwangsaensis and the sister species H. larutensis + H. tehtarik. The ML phylogram places H. cf. titiwangsaensis as the poorly supported (UFB: 78) sister species to H. titiwangsaensis complex and H. larutensis + H. tehtarik being the sister species to them. On the other hand, the BI phylogram places H. larutensis + H. tehtarik as the poorly supported (BPP: 0.72) sister lineage to H. titiwangsaensis complex and H. cf. titiwangsaensis is sister to them (Fig. 2). Both analyses recover populations from Langkawi Island and Gunung Jerai as sister populations, with an uncorrected pairwise sequence divergence of 3.6%. Hemiphyllodactylus cicak from Penang Hill is sister to them, with an uncorrected pairwise sequence divergences of 3.9–4.8% and 4.5–5.7%, respectively (Table 4).

In the nuDNA phylogram, both ML and BI recovered similar topologies (Fig. 3) such as the reciprocally monophyletic typus and harterti groups (1.00/99), as well as the two major clades within the harterti group (0.99/92)—similar to the concatenated phylogram with strong nodal supports. However, the relationships within each clade of the harterti group are poorly resolved and are essentially polytomies.

The BEAST phylogram has an identical topology to the BI phylogram (Fig. 4). All major nodes in the harterti group are generally strongly supported with BPPs of 1.00, except for H. harterti + (H. cicak + (Langkawi population + Gunung Jerai population)) and (H. larutensis + H. tehtarik) + H. titiwangsaensis complex, which are poorly supported with BPP of 0.54 and 0.78, respectively. The divergence times recovered have similar mean ages but slightly broader 95% highest posterior density (95% HPD) intervals than those estimated in Grismer et al. (2015a) and Agarwal et al. (2019). As two Bayesian analyses (BI and BEAST) were consistent in recovering similar topologies it was preferred over the ML topology for ancestral range reconstruction as opposed to forcing the ML topology in a separate BEAST analysis. A separate BEAST analysis that forced the monophyly of H. titiwangsaensis and H. cf. titiwangsaensis> was also performed, however, it produced incorrect topologies for several other species. Therefore, the BEAST analysis that did not force monophyly was considered as more reliable. In the BioGeoBEARS analysis, the DIVALIKE+J model is selected as the best-fit model with the lowest AIC and highest likelihood scores (Fig. 5; Table S1).

The BioGeoBEARS analysis indicates that the ancestor of the harterti group evolved in Peninsular Malaysia approximately 44 (59–34) million years ago (Mya) and diverged into two distinct populations approximately 29 (43–18) Mya, one in the central region in the Banjaran Titiwangsa, and the other most probably in the east, at the Banjaran Timur. At approximately 18 (28–11) Mya, the population in the Banjaran Titiwangsa diverged to form an ancestral population of (H. larutensis + H. tehtarik) + H. titiwangsaensis complex, and another population that diverged into H. cf. titiwangsaensis —both of which remained in the Banjaran Titiwangsa. The ancestral population of (H. larutensis + H. tehtarik) + H. titiwangsaensis complex subsequently diverged at approximately 16 (25–10) Mya, giving rise to H. titiwangsaensis complex in the Banjaran Titiwangsa and another population that diverged later approximately 5 (10–2) Mya into H. larutensis in the west, in the Banjaran Bintang, and H. tehtarik in the east, in the Banjaran Timur. The ancestor of this population likely occurred somewhere between their current distributions on opposite sides of the peninsula.

The other population in the Banjaran Timur diverged approximately 11 (18–7) Mya into H. bintik in the east, and an ancestral population of H. harterti and species found along the northwestern coast. Similarly, the ancestor of this population likely occurred somewhere between their current distributions on opposite sides of the peninsula. The ancestor of H. harterti and species found along the northwestern coast subsequently diverged approximately 9 (15–5) Mya to form H. harterti in the Banjaran Bintang and the ancestor of island species and Gunung Jerai population from the northwestern coast. At approximately 3 (6–2) Mya, the ancestor of the northwestern coast populations diverged to form H. cicak on Penang Island, and the sister populations on Langkawi Island and Gunung Jerai diverging at approximately 2 (4–1) Mya.

The morphological data taken from the Hemiphyllodactylus specimens examined in this study are available in Tables S2–S5. The comparisons of the characters across all species in the harterti group are in Table 5. Populations from Langkawi Island and Gunung Jerai showed distinct morphological differences from all other species of the harterti group (Table 5). The results of the ANOVA, Welch’s ANOVA, and Kruskal-Wallis tests for 10 meristic characters and six size-corrected morphometric characters, as well as the respective post hoc tests for characters with significant mean differences are available in Tables S6 and S7. Based on the results, Langkawi population has significantly higher mean values for SL, IL, FL1, TL1, SVL, AG, HL, SN, HW and ED, and a lower mean value for VS from H. cicak. From the Gunung Jerai population, the Langkawi population shows significantly higher mean values for DS, FL1, TL1, and a significantly lower mean value for AG. Gunung Jerai specimen shows significantly higher values for AG, and HW, and significantly lower values of DS, and VS from H. cicak. The range, mean, and degree of differences among these characters are presented in Table 6 and visualised in Figure 6.

Based on the results of MFA, Langkawi population shows distinct separation from Gunung Jerai population and H. cicak (Fig. 7A). For Dim 1, morphometric characters contributed to 50.6% of the variation, followed by meristic characters (Fig. 7B). For Dim 2, meristic characters contributed to 84.6% of the variation, followed by morphometric characters (Fig. 7B). See Table S9 for more details on the MFA results. The PERMANOVA, demonstrated that the centroid positions between Langkawi population and H. cicak are significantly different (p value = 0.03), but that between Gunung Jerai population with both Langkawi population and H. cicak are not significantly different (p value = 0.25 and 0.2 respectively) likely due to an n = 1 for the former (Table S10).

Based on the results of the phylogenetic and morphological analyses, the populations from Langkawi Island and the Gunung Jerai specimen are described as new species below.

The sky-islands of Peninsular Malaysia harbour high biodiversity and many endemic species have been described (e.g., Grismer et al. 2014c, 2016a, 2017; Davis et al. 2016; Quah et al. 2017a, 2017b, 2020). Some of these species were once thought to be a single, widespread species, which were later found to be composed of multiple species that are endemic to specific localities (e.g., Grismer et al. 2012, 2014c). Besides that, there are also some new species which were described based on completely new discoveries during field surveys, and even from some of the most well-studied localities (e.g., Grismer et al. 2015b, 2016a, 2017; Davis et al. 2016; Quah et al. 2020). These discoveries have indicated that many mountain tops in Peninsular Malaysia have their own endemic species. For example, in addition to Hemiphyllodactylus, different species from amphibians (Ansonia Stoliczka, 1870), gekkonids (Cnemaspis, Cyrtodactylus Gray, 1827), agamids (Pseudocalotes Fitzinger, 1842), skinks (Larutia Böhme, 1981, Tytthoscincus Linkem, Diesmos & Brown, 2011), and snakes (Macrocalamus Günther, 1864) can be found on different mountain tops across Peninsular Malaysia (Grismer et al. 2011, 2012, 2014c, 2016a, 2017; Davis et al. 2016; Quah et al. 2020). Some of the congeneric species are even found in sympatry, where they exploit different microhabitats to avoid competition (Grismer et al. 2014c, 2015a, 2017; Quah et al. 2020). Therefore, the discovery of new populations of Hemiphyllodactylus from additional mountain tops of Peninsular Malaysia in this study is not unexpected, and more new populations are likely to be found in the future.

The diversification of the harterti group across the sky-islands of Peninsular Malaysia is clearly not a recent event but a result of climatic fluctuations occurring since the Early Cenozoic. A similar pattern is also seen in taxa with limited dispersal abilities such as the upland anuran species of Ansonia which diversified across Peninsular Malaysia since the Middle Miocene (Grismer et al. 2016b). The rapidly changing environmental conditions during the Pleistocene could be challenging for them, hence, remaining isolated on respective uplands. In contrast, other taxa with high dispersal abilities such as the upland snake species Asthenodipsas lasgalenensisLoredo et al., 2013, A. vertebralis (Boulenger, 1900) and Lycodon butleri Boulenger, 1900, are found on several isolated mountain tops, and with low genetic divergences (<1.0%) among different populations (Loredo et al. 2013; Grismer et al. 2014b). This is suggested to be the result of the rapid, cyclical climatic fluctuations that occurred during the Pleistocene, where these species migrated downslope together with the expanding montane forest during glacial periods, allowing gene flow among different populations, and then became isolated again during interglacial periods. The high dispersal ability of these species resulted in multiple secondary contact throughout the rapid, cyclical climatic fluctuations of the Pleistocene (~50 cycles). As a result, their last isolation is very recent, after the LGM ended around 11,000 years ago, when the Earth entered current interglacial period (Woodruff 2010; Loredo et al. 2013; Grismer et al. 2014b).

Previous study by Cobos et al. (2016) showed that H. titiwangsaensis is a species complex where the northern population from Cameron Highlands and the southern populations from Fraser’s Hill and Genting Highlands have high genetic divergence, but no morphological differences were found between them. In this study, nuDNARAG1 was used to test the hypotheses by Cobos et al. (2016). However, the nuDNA fragment was insufficient to resolve phylogenetic relationships (Fig. 3). Analysis using restriction site-associated DNA sequencing (RADseq) data is ongoing in order to resolve this issue to check for gene flow. Besides that, while examining the relationship within H. titiwangsaensis complex, another new population (H. cf. titiwangsaensis) from Genting Highlands in the Banjaran Titiwangsa was recovered (Fig. 2). It was initially thought of as being another population of H. titiwangsaensis based on morphological appearances. However, there is a notable 14.0–16.5% uncorrected pairwise sequence divergence between this population and H. titiwangsaensis (Table 4). Besides that, the two Bayesian analyses (BI and BEAST) recovered that this population is not even the closest relative of H. titiwangsaensis, while ML analysis recovered this population as sister to H. titiwangsaensis (Figs 2, 4). The results of BioGeoBEARS suggest that this population might be the intermediate species between H. larutensis in the west and H. tehtarik in the east (Fig. 5). Additional analyses and more detailed examination are required to investigate the relationship between this population and H. titiwangsaensis, and such work is currently in progress.

Within the harterti group, two species occur in sympatry in the Banjaran Bintang (H. harterti and H. larutensis), and another two species occurs in sympatry in the Banjaran Timur (H. bintik and H. tehtarik) (Fig. 1). However, the two sympatric species are not each other’s closest relatives, instead, H. harterti is more closely related to H. bintik, and H. larutensis is more closely related to H. tehtarik (Fig. 2). As eastern and western Peninsular Malaysia are separated by the Banjaran Titiwangsa, it is highly possible that there are yet-to-be discovered populations present in the Banjaran Titiwangsa. These populations would act as the intermediate species between H. bintik found in the east and H. harterti + northwestern coast species found in the west. This would match the topology and geography of the other clade where an ancestral population in the Banjaran Titiwangsa gave rise to H. tehtarik in the east and H. larutensis in the west.

The climate of Sundaland during the Early Cenozoic was generally warm and wet, with the Late Palaeocene Thermal Maximum and the Early Eocene Climatic Optimum occurring around 56 Mya and 52 to 50 Mya, respectively (Zachos et al. 2001; Morley 2012). The global climate started to cool around the Middle Eocene and eventually entered a global icehouse condition around the Early Oligocene (~34 Mya, the Oligocene isotope Event 1 (Oi-1) glaciation) (Zachos et al. 2001). In Sundaland, the climate during the Early to Middle Oligocene was generally cool and dry, and became warmer and wetter during the Late Oligocene, but with some short intervals of cooler periods (Morley 2012). The warm and wet climate reflects the Late Oligocene warming, which occurred approximately 28 to 25 Mya (Zachos et al. 2001). As the ancestral population of the harterti group evolved somewhere in Peninsular Malaysia around 44 Mya, these climatic events might have forced them to move to cooler upland refuges, resulting in their isolation and divergence into two different populations approximately 29 Mya, one in the Banjaran Titiwangsa, and another in the Banjaran Timur. During the Early to Middle Miocene, the climate was generally warmer and wetter, but there was a short period of cooler climate around 23 Mya, during the boundary of Oligocene and Miocene (the Miocene isotope Event 1 (Mi-1) glaciation) (Zachos et al. 2001; Morley 2012). Cooler periods allowed montane forest to expand downslope and the two isolated population might have expanded together with the montane forest to lower altitudes. Then, during the warmer period, montane forest contracted and migrated upslope, creating fragmented montane forests in different upland regions. This might have caused the isolation of the ancestral population of H. cf. titiwangsaensis from the ancestral population of (H. tehtarik + H. larutensis) + H. titiwangsaensis complex at approximately 18 Mya. The latter then gave rise to the ancestral population of H. tehtarik + H. larutensis in the Banjaran Timur and the ancestral population of H. titiwangsaensis complex in the Banjaran Titiwangsa, which diverged around 16 Mya. Although the BioGeoBEARS analysis shows that the ancestral population of H. tehtarik + H. larutensis likely occurred in the Banjaran Timur, it might also be possible that this ancestral population was isolated and diverged from the ancestral population of H. titiwangsaensis complex somewhere within the Banjaran Titiwangsa, then dispersed to the surrounding area during cooler period, and became isolated again on the Banjaran Bintang and the Banjaran Timur, which are separated by the Banjaran Titiwangsa.

The warm and wet climate in Sundaland peaked during the Middle Miocene, when the Mid-Miocene Climatic Optimum occurred around 17 to 15 Mya (Zachos et al. 2001; Morley 2012). According to Morley (2012), intervals of cooler climate occurred in Sundaland, starting from the Latest Middle Miocene, indicating the initiation of the Late Neogene cooling. Since then, the climate fluctuated more frequently between warmer and cooler period, with at least 13 fluctuating cycles occurring during the Middle Late Miocene (9 to 7 Mya) and about six occurring during the Later Late Miocene to Early Pliocene. The climatic fluctuations intensified from the Late Pliocene and continued through the Pleistocene (Morley 2012). The sequential evolution of H. titiwangsaensis complex might reflect fluctuating climates during the Middle Miocene to Pliocene, where montane forests expanded downslope during cooler periods and became continuous throughout the Banjaran Titiwangsa. This facilitated dispersal of the ancestral populations along the region, and then became isolated again at different peaks when the climate warmed. The evolution of the clade containing H. harterti, H. bintik and the northwestern coast species and the divergence between H. larutensis and H. tehtarik might also be due to the climatic fluctuations that occurred between the Oligocene and Pleistocene. Fluctuating sea level associated with climatic fluctuations during the Miocene and Pliocene likely played an important role in the evolution of the northwestern coast species. The water depth around Langkawi and Penang islands is approximately 5 m near the mainland and increased gradually to about 50 m at western side of the islands (Ismail et al. 2012). Global sea level changes during the Oligocene to Early Pliocene were between 30 to 60 m (Miller et al. 2005). During intervals of cooler periods, the sea level became lower, exposing the submerged land between mainland Peninsular Malaysia and the offshore islands. This, together with the downslope expansion of montane forest during the cooler periods, might have led to the dispersal of the ancestral populations between mainland Peninsular Malaysia, Langkawi Island and Penang Island. When the climate turned warm and sea level rose, these populations became isolated, giving rise to H. cicak from Penang Hill, H. puncaksp. nov. from Langkawi Island, and H. jeraiensissp. nov. from Gunung Jerai.