*Article* **Karstic Landscapes Are Foci of Species Diversity in the World's Third-Largest Vertebrate Genus** *Cyrtodactylus* **Gray, 1827 (Reptilia: Squamata; Gekkonidae)**

**Lee Grismer 1,\*, Perry L. Wood, Jr. 2, Nikolay A. Poyarkov 3,4, Minh D. Le 5,6,7, Suranjan Karunarathna 8, Siriwadee Chomdej 9, Chatmongkon Suwannapoom 10, Shuo Qi 11, Shuo Liu 12, Jing Che 13, Evan S. H. Quah 1,14, Fred Kraus 15, Paul M. Oliver 16, Awal Riyanto 17, Olivier S. G. Pauwels <sup>18</sup> and Jesse L. Grismer <sup>1</sup>**


**Abstract:** Karstic landscapes are immense reservoirs of biodiversity and range-restricted endemism. Nowhere is this more evident than in the world's third-largest vertebrate genus *Cyrtodactylus* (Gekkonidae) which contains well over 300 species. A stochastic character mapping analysis of 10 different habitat preferences across a phylogeny containing 344 described and undescribed species recovered a karst habitat preference occurring in 25.0% of the species, whereas that of the other eight specific habitat preferences occurred in only 0.2–11.0% of the species. The tenth category—general habitat preference—occurred in 38.7% of the species and was the ancestral habitat preference for *Cyrtodactylus* and the ultimate origin of all other habitat preferences. This study echoes the results

**Citation:** Grismer, L.; Wood, P.L., Jr.; Poyarkov, N.A.; Le, M.D.; Karunarathna, S.; Chomdej, S.; Suwannapoom, C.; Qi, S.; Liu, S.; Che, J.; et al. Karstic Landscapes Are Foci of Species Diversity in the World's Third-Largest Vertebrate Genus *Cyrtodactylus* Gray, 1827 (Reptilia: Squamata; Gekkonidae). *Diversity* **2021**, *13*, 183. https://doi.org/ 10.3390/d13050183

Academic Editor: Eric Buffetaut

Received: 23 March 2021 Accepted: 23 April 2021 Published: 28 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

5

of a previous study illustrating that karstic landscapes are generators of species diversity within *Cyrtodactylus* and not simply "imperiled arks of biodiversity" serving as refugia for relics. Unfortunately, the immense financial returns of mineral extraction to developing nations largely outweighs concerns for biodiversity conservation, leaving approximately 99% of karstic landscapes with no legal protection. This study continues to underscore the urgent need for their appropriate management and conservation. Additionally, this analysis supports the monophyly of the recently proposed 31 species groups and adds one additional species group.

**Keywords:** Indochina; Southeast Asia; phylogeny; Indo-Australian Archipelago; Bent-toed geckos; karst; conservation

#### **1. Introduction**

The dramatic topography of karstic landscapes composes some of the most surreal images of our world and has stirred the emotions of ancient artisans and natural historians for time on end. But not only are these crenulated, repeating layers of rugged terrain steeped in natural beauty (Figure 1), they are the only refuge for some of the most seriously endangered species on the planet [1]. Asia contains 8.35 million km2 of karstic habitat with some of the most extensive concentrations ranging from China to western Melanesia (Figure 2). These formations are notable for their fragmented, island-like nature, with hills, caves, and towers forming archipelagos of habitat-islands stretching across broad geographic areas. This, and their fractured and eroded surfaces—which provide a myriad of microhabitats in which many taxonomic groups have specialized—have contributed to their extraordinarily high degrees of range-restricted endemism [2–5]. Karst formations are often referred to as "imperiled arks of biodiversity" [5]. However, a stochastic character mapping analysis of habitat preference using 243 species of the gekkonid genus *Cyrtodactylus* the third most speciose vertebrate genus on the planet—indicated just the opposite [6]. Grismer et al. [6] demonstrated that karstic landscapes not only harbor range-restricted endemics, but have been the foci of speciation for the largest independent gekkonid radiations across all of Indochina and Southeast Asia. They went on to show that even in this ecologically labile genus, karst-associated species outnumbered by threefold all other species bearing other specific habitat associations. As such, this has transformed our view of karstic landscapes from that of "limestone museums" harboring relictual endemics, to platforms of speciation and generators of biodiversity across a broad *taxonomic* landscape (e.g., [7–10]).

**Figure 1.** The Dragon Back karst formation in Perlis State, Peninsular Malaysia.

**Figure 2.** The distribution of karstic landscapes throughout Indochina and the Indo-Australian Archipelago.

*Cyrtodactylus* is by far the most speciose and ecologically diverse gekkotan genus [6,11]. It currently contains 306 nominal species (as of 14 February 2021; [12]) ranging from South Asia to Melanesia (Figure 3) where they occupy a vast diversity of habitats. As would be expected from a group this large and widely distributed, it bears a broad variety of ecotypes ranging from short robust terrestrial species to cryptically colored arboreal species to gracile cave-dwelling and karst-adapted specialists (e.g., [13–21]; Figure 4). The annual rate at which new species are being described is unprecedented and shows no signs of leveling off (Figure 5) and the majority of the most recently described species are associated with karst formations. In some cases, multiple species from distantly related clades may be found throughout a single karstic archipelago [14,20], and even more remarkable, different species from distantly related clades may even occupy the same small karst formation [14,20]. The intent of this paper is to test, (1) whether or not the same clades bearing the same specific habitat preferences presented by Grismer et al. [6] are recoverable, (2) whether or not the relative frequencies of species in each habitat preference category are not significantly different than that reported by Grismer et al. [6], (3) and specifically, is the hypothesis that karstic landscapes are generators of biodiversity further supported. We test these hypotheses by augmenting Grismer et al.'s [6] original phylogeny of 243 species with an additional 101 species (a 44% increase in species coverage) and by adding a new category of habitat preference. Additionally, with this significant influx of species, we test the monophyly of the 31 different species groups recently designated by Grismer et al. [11] based on their phylogeny of 310 named and unnamed species (an 11% increase in species coverage).

**Figure 3.** Generalized distribution of the genus *Cyrtodactylus*.

**Figure 4.** Representative ecotypes of the 10 different habitat preferences in the genus *Cyrtodactylus.* Photographs by (**A**) L. Lee Grismer, (**B**) Steve J. Richards, (**C**,**D**) L. Lee Grismer, (**E**) Suranjan Karunarathna, (**F**) L. Lee Grismer. (**G**) Evan S. H. Quah, (**H**,**I**) L. Lee Grismer, and (**J**) Peter Geissler.

**Figure 5.** Cumulative number of species of *Cyrtodactylus* described per year. The trajectory of new species descriptions from 2000 to 18 March 2021 indicates that the true diversity of this genus is not yet calculable and that 48% of the newly described species during this period have come from Myanmar, Vietnam, and Malaysia.

#### **2. Materials and Methods**

#### *2.1. Habitat Preferences and Ecotypes*

Here we refine some of the criteria for designating habitat preference used by Grismer et al. [6] based on newly acquired data from recent publications and fieldwork. We also add an additional habitat preference (sandstone), bringing the total to 10 as opposed to nine categories (Table S1). Habitat preference for each species was coded as a discrete character state and ascertained by integrating data from the literature, firsthand experience of the authors, and personal communication with researchers familiar with particular species. Grismer et al. [6] acknowledged that some of these categories could be further subdivided (e.g., arboreal into branch, twig, and leaf), but those subdivisions become far less defensible owing to a lack of detailed microhabitat information. In this regard, many species can be considered data deficient, inasmuch as baseline information on their ecological requirements are often limited to anecdotal observations made at the time of their collection (e.g., [14]). The potential biases of using limited observations from a single locality at one point in time to ascertain the habitat preference of an entire species does not go unnoticed. However, in many cases, these are the only data available. Nonetheless, judiciously vetted, natural history observations summarized across the literature coupled with our own field observations and those of others, can provide a useful framework for supporting robust, testable, downstream hypotheses regarding habitat preference. The habitat preferences and their associated ecotypes bearing the same categorical names are described below. Obvious morphological correlates associated with some ecotypes are noted only for additional clarity.

1. General (Figure 4A). Species that use the majority of the microhabitats in their immediate surroundings in whatever environment they inhabit. The microhabitats may include rocks of all types (when present), logs, tree trunks (with or without holes and crevices), and all vegetative structures of various dimensions, the ground, and humanmade structures in many cases. No particular microhabitat is notably preferred over any other although some species may be most often observed in low vegetation.


This species is similar in body shape to closely related granite-associated species (Grismer unpublished).

#### *2.2. Mitochondrial DNA*

The data set of Grismer et al. [6] was augmented with 107 additional ingroup species resulting in a matrix composed of 344 described and undescribed species (i.e., species identified in previous phylogenies but not yet described) of *Cyrtodactylus* (Table S1). A phylogeny was constructed using 1474 base pairs of the mitochondrial gene NADH dehydrogenase subunit 2 gene and its flanking tRNAs (hereafter referred to as ND2). *Agamura persica*, *Bunopus tuberculatus*, *Hemidactylus angulatus*, *H. frenatus*, *H. garnotii*, *H. mabouia*, *H. turcicus*, *Lialis jicari*, *Mediodactylus russowii*, *Mokopiriakau cryptozoicus*, *Pygopus nigriceps*, *Sphaerodactylus torrei*, *Stenodactylus petrii*, *Tenuidactylus elongatus*, *Toropuku stephensi*, and *Tropiocolotes steudneri*—encompassing all other major gekkotan lineages—were used to root the tree following Wood et al. [38]. Genomic DNA was isolated from liver or skeletal muscle from new tissue samples stored in 95% ethanol, using standard phenol-chloroformproteinase K (final concentration 1 mg/mL) extraction procedures with subsequent isopropanol precipitation following Hillis et al. [39] or a SPRI magnetic-bead extraction protocol (https://github.com/phyletica/lab-protocols/blob/master/extraction-spri.md; accessed on 15 January 2021). The ND2 gene, with parts of adjacent tRNAs, was amplified using a double-stranded Polymerase Chain Reaction (PCR) under the following conditions: 1.0 μL genomic DNA (10–30 μg), 1.0 μL light-strand primer (concentration 10 μM), 1.0 μL (H5934, 5- – AGRGTGCCAATGTCTTTGTGRTT–3- , following [6]), heavy-strand primer (concentration 10 μM), (L4437b, 5- –AAGCAGTTGGGCCCATRCC–3- , following [6]) 1.0 μL dinucleotide pairs (1.5 μM), 2.0 μL 5 buffer (1.5 μM), MgCl 10× buffer (1.5 μM), 0.1 μL Taq polymerase (5 u/μL), and 6.4 μL ultra-pure H2O. PCR reactions were executed on Bio-Rad T100™ gradient thermocycler 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 55 ◦C for 35 s, followed by a cycle extension at 72 ◦C for 35 s, for 31 cycles. All PCR products were visualized using 1.0% agarose gel electrophoresis. Successful PCR products were sent to Evrogen® (Moscow, Russia), Genetech Sri Lanka Pvt. Ltd. (Colombo, Sri Lanka), or Genewiz® (South Plainfield, NJ, USA) for PCR purification, cycle sequencing, sequencing purification, and sequencing using the same primers as in the amplification step. Sequences were analyzed from both the 3 and the 5 ends separately to confirm congruence between reads. Forward and reverse sequences were uploaded and edited in GeneiousTM 2019.0.4 (https://www.geneious.com). Following sequence editing, the protein-coding region and the flanking tRNAs were aligned using the MAFTT v7.017 [40] plugin under the default settings in Geneious™ 2019.0.4 (https://www.geneious.com). Mesquite v3.04 [41] was used to calculate the correct amino-acid reading frame and to confirm the lack of premature stop codons in the ND2 portion of the DNA fragment.

#### *2.3. Phylogenetic Analyses*

A Maximum likelihood (ML) analysis was implemented using the IQ-TREE webserver [42,43] preceded by the selection of substitution models using the Bayesian Information Criterion (BIC) in ModelFinder [44] which selected TVM+F+I+G4 for the tRNAs and codon position 1 and GTR+F+I+G4 for codon positions 2 and 3. One-thousand bootstrap pseudoreplicates via the ultrafast bootstrap (UFB; [45]) approximation algorithm were employed, and nodes having UFB values of 95 and above were considered strongly supported [46]. We considered nodes with values of 90–94 as well supported.

A Bayesian inference (BI) phylogeny was estimated using Bayesian Evolutionary Analysis by Sampling Trees (BEAST) version 2.4.6 [47] implemented in CIPRES (Cyberinfrastructure for Phylogenetic Research; [48]). Input files were constructed in Bayesian Evolutionary Analysis Utility (BEAUti) version 2.4.6 using a lognormal relaxed clock with unlinked site models, linked trees and clock models, and a Yule prior and run in BEAST version 2.4.6 [47] on CIPRES. bModelTest, implemented in BEAST, was used to

numerically integrate over the uncertainty of substitution models while simultaneously estimating phylogeny using Markov chain Monte Carlo (MCMC). MCMC chains were run for 350,000,000 generations and logged every 35,000 generations. The BEAST log file was visualized in Tracer v. 1.6.0 [49] to ensure effective sample sizes (ESS) were well above 200 for all parameters. A maximum clade credibility tree using mean heights at the nodes was generated using TreeAnnotator v. 1.8.0 [50] with a burn-in of 1000 trees (10%). Nodes with Bayesian posterior probabilities (BPP) of 0.95 and above were considered strongly supported [51,52]. We considered nodes with values of 0.90–0.94 as well supported.

Grismer et al. [6] demonstrated that in their 243-species data set, the third codon position contributed significantly to the strongly supported topological resolution of the tree and showed no signs of codon saturation. In their 310-species tree, Grismer et al. [11] demonstrated that their mito-nuclear tree constructed from ND2 and three nuclear genes did not improve the resolution or the nodal support of the deep nodes in their ND2 tree. Therefore, only ND2 was used in this analysis.

#### *2.4. Ancestral State Reconstruction*

In order to estimate the probability of each habitat preference at each node in the tree, we employed a stochastic character mapping (SCM) analysis implemented in R [v3.4.3] using the R package Phytools [53] on the BEAST tree converted to newick format. The transition-rate matrix that best fit the data was identified by comparing corrected Akaike Information Criterion (AICc) scores among alternate models using the R package ape 5.2 [54]. Three transition-rate models were considered: a 90-parameter model having different rates for every transition type (the ARD model); a 45-parameter model with equal forward and reverse rates between states (the symmetrical rates SYM model); and a single-rate parameter model that assumes equal rates among all transitions (ER). Lastly, an MCMC approach was used to sample the most probable 1000-character histories from the posterior using *make.simmap*() and then summarized them using the *summary*() command.

#### **3. Results**

The ML analysis recovered essentially the same well to strongly supported tree (Figure 6) recovered in Grismer et al. [11]. The same 31 monophyletic species groups designated in Grismer et al. [11] were recovered here even though sampling in was greatly expanded with additional species (Table S1). The ML analysis also recovered a new clade, designated here as the *tibetanus* group, that is composed of *Cyrtodactylus tibetanus*, *C.* cf. *tibetanus,* and *C. zhaoermii*. *Cyrtodacylus* cf. *tibetanus* and *C. zhaoermii* were unavailable for the analysis of Grismer et al. [11], where *C. tibetanus* was recovered as the earliest diverging member of the *lawderanus* group. Che et al. [55] recovered the same new clade in a less inclusive (i.e., fewer species) mito-nuclear phylogeny. Although Grismer et al. [11] recovered *C. rubidus* as the sister species of the *lateralis* group, it was not included in that group because this relationship was well supported only in the ML analysis and not the BI analysis. Here, it is placed in the *lateralis* group with high support in both analyses (90 UFB, 0.90 BPP), a grouping also supported by the fact that all members of this group have prehensile tails.

**Figure 6.** Majority-rule consensus tree from ML bootstrap replicates of 344 species of Cyrtodactylus. Phylogeny based on 1474 base pairs of the mitochondrial gene ND2 illustrating the designation of 32 monophyletic species groups.

The BEAST analysis recovered a tree with generally strong nodal support throughout with a 94.4% topological consistency (recovering 322 of the same 347 nodes) as the ML tree (Figure 7). The AICc scores for the three transition-rate models were ARD = 1101.751; SYM = 1035.445; and ER = 890.9552. The results of the SCM analysis were consistent with

those of Grismer et al. [6] in that the ER model recovered large and small clades that independently evolved the same habitat preferences throughout the geographic range of the genus (Figure 7A). The SCM recovered a general habitat preference as being ancestral for not only the genus *Cyrtodactylus* but for all other major clades and ultimately all other habitat preferences as well. Notably for this study, however, the two largest independently evolved lineages of karst-associated species—the lineage composed of the *sadansinensis*, *yathepyanensis*, *oldhami*, *sinyineensis*, and *chauquangensis* groups and the *angularis* group—were also recovered, even with their expanded species contents. Their parapatric distributions across much of western and northern Indochina coincide with regions bearing the most extensive karstic landscapes (Figure 8). Other less diverse, independently evolved karstic lineages, such as the *linnwayensis* group from the Shan Plateau in Myanmar and a karst-associated subclade from the Thai-Malay Peninsula in the *pulchellus* group, were also recovered and associated with regions rich in karstic habitats (Figures 7A and 8). Several isolated instances of the independent evolution of karst habitat preference are scattered across the tips of the tree, representing species from Borneo (*C. cavernicolus*, *C. limajalur*, *C. muluensis*), Cambodia (*C. laangensis*), China (*Cyrtodactylus* sp. SYS r1232), Indonesia (*C. darmandvillei*), Myanmar (*C. aunglini*, *C. chrysopylos*, *C. myaleiktaung*), Papua New Guinea (*C. tanim*), Peninsular Malaysia (*C. evanquahi*, *C. guakanthanensis*, *C. gunungsenyumensis*, *C. metropolis*, *C. lenggongensis*, *C. sharkari*), and Vietnam (*C.* sp. nov., *C. yangbayensis*) (Figure 7A).

Node probability of ancestral habitat preference Percentage of species in each habitat

**Figure 7.** (**A**) Stochastic character map of the 10 habitat preferences on a maximum clade credibility BEAST phylogeny showing the probability of the ancestral habitat preference at each node and the major clades of karst-associated species groups and their general geographic distribution. (**B**) The habitat preference of each species, and the number of species in each habitat preference category. (**C**) Pie chart showing the percentage of species bearing each of the 10 habitat preferences.

**Figure 8.** Distribution of the major clades of the karst-associated species groups throughout Indochina. Inset illustrates their co-distributions, with the geographic areas bearing the most extensive karstic landscapes.

These data are consistent with those of Grismer et al. [6] in showing that the frequency of karst-associated species far out-numbers that of any other specific habitat preference and is nearly two and one-half times more prevalent than any other specific habitat preference in that it contains 25.0% of the species followed by trunk (11.0%), granite (9.2%), terrestrial (8.4%), arboreal (3.8%), cave (2.0%), swamp (1.4%), and intertidal and sandstone (0.2%; Figure 7B). In Grismer et al. [6], granite-associated species comprised the second highest habitat preference and trunk-associated species the third. That ranking has been reversed here. The percentage of species with a karst habitat preference was 29.6% in Grismer et al. [6] but dropped to 25.0% here. We posit that this drop of nearly 5% is a direct result of our inability to explore unsurveyed karstic regions on the Shan Plateau and in the Salween Basin of Myanmar during 2020 due to COVID-19.

#### **4. Discussion**

The analysis presented here is based on the most complete phylogeny of the genus *Cyrtodactylus* to date with an increase of 101 species from that of Grismer et al. [6] and 35 from that of Grismer et al. [11]. The hypotheses marshaled by Grismer et al. [6] concerning the evolution of habitat preference is supported here in that there was no notable change in the frequencies of species bearing different habitat preferences across the genus—even with the addition of 107 species. More specifically, however, a karst habitat preference retained a higher frequency than that of any other specific habitat preference (25.0% versus 0.2–11%), supporting the hypothesis that these landscapes are platforms for the generation of biodiversity. This pattern is particularly strong in Indochina and less so on islands throughout the Indo-Australian Archipelago, reflecting the sharp contrast in the extent of karstic landscapes between these regions (Figure 8). These data clearly underscore the importance of karstic habitats to this hyper-diverse genus and continue to amplify the work of many other authors indicating that the high levels of biodiversity and rangerestricted endemism in karstic habitats rivals that of most other habitats throughout the tropics (see discussions in [1,4,5,10,56–61]). The sad irony is that, although these are some of the most imperiled ecosystems on the planet due to unregulated and unsustainable quarrying practices, only 1% of these terrains throughout Asia are afforded any form of

legal protection. Therefore, the diversity of the karst-associated species in general—and *Cyrtodactylus* in particular—are, for the most part, without legal protection. Unfortunately, the immense financial returns from cement manufacturing makes the challenge of karst conservation difficult and many governments from developing nations that are willing to overlook sustainable quarrying policies in order to expand their economy [1]. Continued exploitation of karstic habitats for limestone shows no signs of abating.

#### **5. Conclusions**

This study echoes the results of Grismer et al. [6] in that karstic landscapes are exceedingly important for maintaining *Cyrtodactylus* diversity and serve as foci for their speciation and maintenance of their diversity. Referring to them as "imperiled arks of biodiversity" is somewhat misleading as these are ecological platforms for speciation that not only continue to generate the most speciose, independent, radiations of the Gekkota, but do so across a broad range of other taxonomic groups (e.g., [7–10,62]). Referring to them as "imperiled arks of biodiversity" instead of centers for speciation draws attention away from their importance as generators of biodiversity in an era of biodiversity crisis and could potentially lessen the urgency for legislative conservation measures.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/d13050183/s1, Table S1: Species, habitat preference with supporting references, species group designations, and GenBank accession numbers for specimens used in the SCM analysis. Species can be cross-referenced to Figure 6 by their GenBank no.

**Author Contributions:** Conceptualization, L.G., J.L.G. and P.L.W.J.; methodology, L.G., J.L.G. and P.L.W.J.; formal analysis, L.G. and P.L.W.J.; resources, N.A.P., M.D.L., S.K., S.C., C.S., S.Q., S.L. and J.C.; data curation, L.G. and P.L.W.J.; writing—original draft preparation, L.G.; writing—review and editing, all authors.; visualization, L.G.; supervision, L.G.; funding acquisition, S.K. and N.A.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially supported by the Nagao Natural Environment Foundation (2018- 20) grant to Suranjan Karunarathna for field and lab works. This research was partially funded by the Russian Science Foundation grant number 19-14-00050 to Nikolay A. Poyarkov (sample collection, molecular and phylogenetic analyses). This research was funded by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under the grant number 106.05-2020.24 and the National Geographic grant number 230151 to Minh Duc Le.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** This paper follows all the relevant guidelines and ethics at https: //www.mdpi.com/ethics.

**Acknowledgments:** We thank Steve Richards for additional information and photographs of taxa from New Guinea. N.A. Poyarkov thanks R.A. Nazarov, Valentina F. Orlova, Evgeniya N. Solovyeva and Eduard A. Galoyan (ZMMU, Russia), Andrei N. Kuznetsov and Svetlana P. Kuznetsova (JRVTC, Hanoi, Vietnam), Nikolai L. Orlov and Natalia B. Ananjeva (ZISP RAS, Russia), Parinya Pawangkhanant (AUP, Thailand), Than Zaw and May Thu Chit (Mandalay University, Myanmar), Peter Brakels (IUCN, Laos), Nguyen Van Tan (SVW, Vietnam), Indraneil Das (UNIMAS, Malaysia), Hung Ngoc Nguyen and Si-Min Lin (NTNU, Taiwan), Anna B. Vassilieva (IPEE RAS, Russia), Platon V. Yushchenko, Anna S. Dubrovskaya, Sabira S. Idiatullina, and Vladislav A. Gorin (MSU, Russia) for continuous support and assistance in the field and in the lab. S. Karunarathna thanks Kanishka Ukuwela, Anslem de Silva, Majintha Madawala, Dinesh Gabadage and Madhava Botejue for various support and assistance in the field and in the lab and the Department of Wildlife Conservation (WL/3/2/42/18a, b, c) and Forest Conservation Department (RandE/RES/NFSRCM/2019-04) for research permits. Chatmongkon Suwannapoom thanks the Unit of Excellence 2020 on Biodiversity and Natural Resources Management, University of Phayao (UoE63005). PLW's contribution to this paper constitutes contribution number 943 of the Auburn University Museum of Natural History.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

