Phylogeography and ecological niche modeling reveal evolutionary history of Leiolepis ocellata (Squamata, Leiolepidae)

Abstract Leiolepis ocellata is a lizard species distributing in topographically diverse habitats in northern Thailand. To explore its evolutionary history, 113 samples of L. ocellata were collected from 11 localities covering its distributional range in northern Thailand, and sequenced for mtDNA fragments (Cyt b and ND2). Pairwise comparisons across sampling localities yielded significant genetic differentiation (F ST and Jost's D) but no clear pattern of isolation by distance could be demonstrated based on the Mantel test. Phylogenetic and network analyses highlighted six haplogroups. Their divergence times were estimated to occur during the Pleistocene, much more recent than major orogenic events affecting northern Thailand. Instead, the results suggested that lineage divergences, of particularly eastern and western haplogroups of the region, coincided with the major rivers in the region (Yom river and Ping river, respectively), indicating vicariance in response to riverine barriers. Furthermore, ecological niche modeling suggested an expansion of suitable habitats of L. ocellata, when LGM‐liked conditions. This expansion potentially facilitated their dispersal among adjacent localities leading to lineage diversification and genetic admixture, after the riverine divergence.

. Other environmental factors especially climatic conditions are also important in influencing the dispersal ability of animals, shaping their geographic distributions, and contemporary population structures (Orsini et al., 2013;Pyron & Burbrink, 2010). Incorporating genetic data and species distribution modeling based on paleoclimatic information provides important insights into the historical dispersal patterns of species (Canestrelli et al., 2007;Kim et al., 2020;Moore, 1995;Ujvari et al., 2008;Werneck et al., 2012). The paleoclimatic events during the Pleistocene (2.5 mya-11.7 ka) are interesting since the multiple oscillations of climate, as glacial and interglacial periods, established repeated contraction and expansion of suitable habitat and species distribution range, contributing genetic divergence, and genetic admixture due to a subsequent secondary contact (Canestrelli et al., 2007;Ding et al., 2011;Graham et al., 2013;Grismer et al., 2014;Huang et al., 2013;Lanier & Olson, 2013;Lin et al., 2010;Nicolas et al., 2018).
The Indo-Burma region is an intriguing biological hotspot for phylogeographic studies due to great biodiversity and complex landscapes and climates that are consequences of multiple paleogeographic events (Hall, 2009;Myers et al., 2000). That part in northern Thailand represents an important ecoregion topographically characterized by north-south-oriented mountain ranges and intermontane lowlands.
The mountain ranges are extensions of the Tibetan Plateau, resulting from the accretion of the Indian subcontinent with Eurasia during the Tertiary (~50 mya; Tapponnier et al., 2001;Yin & Harrison, 2000).
The Tibetan plateau is also the origin of large waterways in northern Thailand, including the Salween and Mekong rivers (Brookfield, 1998;Clark et al., 2004;Horton et al., 2002;Nie et al., 2018). Geological and biological evidences indicate that these two rivers formed the Siam river system, flowing via the present-day northern Thai Ping river (the longest river in northern Thailand; Woodruff, 2010), and Yom river (Brookfield, 1998), respectively, into the Gulf of Thailand via the Chao Phraya river valley of central Thailand, during the Pleistocene (Hutchison, 1989;Woodruff, 2010). These historical events not only shaped the geological landscape but also brought a unique climate to northern Thailand promoting habitat heterogeneity (Heaney, 1991;Penny, 2001;White et al., 2004;Woodruff, 2010). These habitats harbor many endemic species and their conservation is of major concern (Bain & Hurley, 2011;Chan-ard et al., 2015;Das, 2010;Trisurat et al., 2015). Despite evolutionary impact of these major historical, geographic, and climatic events, their influence on the phylogeography and population structure of terrestrial fauna confined to northern Thailand has scarcely been investigated.
Leiolepis originated on the Southeast Asian tectonic plate, which separated from the northern margin of the Australia-New Guinea plate and collided with Asia around 120 mya or earlier (Macey et al., 2000).
In this study, we sampled and sequenced L. ocellata across their geographic ranges in northern Thailand and conducted phylogenetic, population genetic analyses using two partial sequences of mitochondrial genes (Cyt b and ND2), and ecological modeling.
The main goals of this study were to explore how paleogeographic events, including orogenic development and river systems, and paleoclimate might explain observed phylogeographic and population structure of L. ocellata populations in northern Thailand.

| Sample collection
In this study, L. ocellata were sampled and identified to species based on morphological characteristics (Arunyavalai, 2003;Peters, 1971) our additional surveys, throughout northern Thailand (Promnun et al., 2020). Additionally, we sampled two individuals of L. belliana from Kamphaeng Phet Province (KP) to use as outgroups. The lizards were randomly captured by hand and traps. Tissues of tail tip (<1 cm in length) were clipped using sterile scissors. Collected tissues were preserved in 95% ethanol and stored at −20°C until further analyses.
Pairwise F ST calculated by Arlequin was used as a genetic distance matrix. Geographic distances among populations were measured from the distance between sampling localities using QGIS (QGIS Development Team, 2020).

| Phylogenetic and network analyses
Phylogenetic trees of L. ocellata collected in this study were con- partitioned allowing the evolutionary model to specify each partition separately. Kakusan4 (Tanabe, 2007) was used to select the best-fit evolutionary model under Akaike information criterion (AIC; Akaike, 1974) and Bayesian information criterion (BIC; Schwarz, 1978) for ML and BI, respectively. The selected models for the ML trees were: Cyt b, J1 Invariant; and ND2, TN93 Gamma; and for the BI trees were: Cyt b, HKY85 Invariant; and ND2, HKY85 Gamma. The ML trees were constructed using RAxML v8.2.12 (Stamatakis, 2014) on CIPRES (Miller et al., 2010) with 1,000 pseudoreplicates to estimate branch confidence values. The bootstrap value of 70% or higher was considered as significant support (Huelsenbeck & Hillis, 1993). The BI trees were constructed in MrBayes v3.2.6 (Huelsenbeck & Ronquist, 2001) under a Metropolis-coupled, Markov chain Monte Carlo (MC-MCMC) approach, started from random tree, run twice in parallel with a fourchain analysis for 10 million generations. The trees were sampled every 100 generations and 25% of the generations were discarded as burnin. The outputs were checked by examining Effective Sample Size (ESS; >200) in Tracer v1.7.1 (Rambaut et al., 2018). The remaining trees were used to estimate the consensus topology, branch length, and posterior TA B L E 1 Sampling localities, sample size (N), number of haplotypes (n), haplotype diversity (h), nucleotide diversity (π) based on concatenated mtDNA probability. The posterior probability of 95% or higher was considered as a significant support (Larget & Simon, 1999). The trees were visualized and edited in FigTree v1.4.3 (Rambaut, 2009). We additionally established Median-Joining Network for mtDNA dataset to illustrate relationships among haplotypes using Popart (Leigh & Bryand, 2015).

| Divergence time estimation
Divergence times were estimated from the Cyt b dataset of L. ocellata and L. belliana. Analyses were conducted in BEAST v1.10.4 (Drummond & Rambaut, 2007) using uncorrelated relaxed clock model with constant-size coalescent prior and the model of nucleotide substitution described above. Since recent fossil calibration for Leiolepis was not available, we used the substitution rates of Cyt b derived from of many species of lizards following previous study of L. reevesii (Lin et al., 2010). We assumed mean substitution rate as a normal distribution, with a mean of 0.01 and a standard devia-

| Ecological niche modeling
Ecological niche modeling was conducted using Maxent v3.4.1 (Phillips et al., 2006) Hijmans et al., 2005). As the preliminary models predicted that the distribution of L. ocellata occurred in northern Thailand, we clipped and only presented the map of northern Thailand using QGIS. To select the most important predictor variables, we ran initial models with the default setting for each period.
We selected the predictor variables from the layer of the present time, which were >10% for percent contribution and permutation importance, and less correlated to each other (Pearson correlation coefficient |r| < 0.80; Khanum et al., 2013). The predictor variables included bio3, bio7, and bio15. The subsequent models were generated with a randomly selected 75% and the remaining 25% of presence data being used as training and testing data, respectively, for model validation (Corbalán et al., 2011). Models were under 10 replications of bootstrap replicated run type (5,000 iterations) and other parameters with default setting. The areas under the curve (AUC) of a receiver operating characteristics (ROC) plot were considered to evaluate performance of models. The higher the AUC values appear, the more reliable the models are (Eskildsen et al., 2013). The contribution of each selected variable was assessed from the percent contribution and permutation importance.

| Population genetic analyses
The final alignments of concatenated partial mitochondrial DNA

| Phylogenetic and network analyses
As concatenated mtDNA trees based on both BI and ML approaches provided similar topology, only the BI tree was shown in Figure 2. The phylogenetic tree topology strongly supported F I G U R E 2 Bayesian Inference phylogenetic trees based on concatenated mtDNA (Cyt b and ND2) of L. ocellata with L. belliana from KP as an outgroup. The two last alphabets indicated code of sampling location following   value and 100% posterior probability) was distributed in central northern Thailand, and found in MP, CS, CM, and TK.
The network analyses of concatenated mitochondrial genes provided similar results to those of the phylogenetic trees but revealed slightly different genealogy (Figure 3). A total 42 haplotypes of L.
ocellata formed six recognizable groups (A-F) identified in the trees.
Each locality harbored private haplotypes and most haplotypes were unique to their geographic localities except for two haplotypes shared between MN and WS.

| Divergence time
Based on Cyt b mtDNA sequences, the coalescence times of the L. ocellata haplogroups were all estimated to have occurred during the Pleistocene (Table 3, Figure 4). The first divergence of L. ocellata (node I in Figure 4) was estimated to at 0.81 mya (HPD = 0.24-1.62).
Divergence within the first group (node II in Figure 4) occurred at 0.61 mya (HPD = 0.14-1.29). Within the second group, divergence between haplogroup C and the remaining haplogroups (node III in

| Population genetics of L. ocellata
Exploring the phylogeography and population genetics of organisms provides valuable insights for both evolutionary processes and conservation of biodiversity (Avise, 2009;Frankham et al., 2010;Manel et al., 2003). The effects of paleogeographic and climatic events on terrestrial animals were widely observed in topographically diverse regions, which often serve as biological hotspots (e.g., Carnaval et al., 2009;Grismer et al., 2014;Guo et al., 2011;Lin et al., 2010;Zhu et al., 2016). Northern Thailand constitutes part of the Indo-Burma biodiversity hotspot, possessing complex geography and high habitat heterogeneity, but has rarely been considered for phylogeographic study. Our investigation on L. ocellata may be the first phylogeographic and population genetic study of any terrestrial fauna confined to northern Thailand.
Evidence from mtDNA enabled us to uncover significant genetic diversity within L. ocellata as indicated by its high haplotype diversity and π = 1.06; Lin et al., 2010). However, the genetic diversity within some populations such as those at DK and MN was relatively low, which could be due to many factors, such as overexploitation as food by local people that caused population decline and reduction in genetic diversity (Arunyavalai, 2003;Lin et al., 2010).
Analyses from the concatenated mtDNA implied strong population genetic structure of L. ocellata populations as indicated by significant genetic differentiation (pairwise Fst and Jost's D) and lack of shared haplotypes among populations other than for the MN-WS. This population structure was likely a consequence of low dispersal ability of lizards (Clark et al., 1999;Dubey & Shine, 2010;Koumoundouros et al., 2009;Orsini et al., 2013). Dispersal of Leiolepis seems to be limited by their behaviors and life history (Arunyavalai, 2003). Based on our field observations, each burrow of the ground-dwelling Leiolepis was usually occupied by one adult (with juveniles in the same burrow for mature females). Juveniles usually dispersed to adjacent areas after maturation and remained in a particular burrow unless they were disturbed. It was also noticed that some adults protected their burrows and adjacent areas against other individuals indicating that territorial behavior that may exist in factors might further attribute to its phylogeographic structure at a regional scale.

| Phylogeographic pattern and divergence between major lineages
Both phylogenetic and network analyses, based on concatenated mtDNA, provided a similar phylogeographic pattern within L. ocellata in which six haplogroups could be identified. To further infer the potential factors related to this phylogeographic pattern, divergence times were estimated from Cyt b tree using substitution rates from other reptiles (Agamidae), following the study of L. reevesii (Lin et al., 2010). Without fossil calibration of Leiolepis, the use of these substitution rates should be interpreted with caution. According to these rates, the divergence times of six major lineages of L. ocellata were all estimated to occur during the Pleistocene. These divergence times were congruent with those of other species, which were suggested to be influenced by several scenarios (e.g., Canestrelli et al., 2007;Gonçalves et al., 2018;Graham et al., 2013;Grismer et al., 2014;Lin et al., 2010;Nicolas et al., 2018). While in general, population genetic structure may be explained by an isolation by distance pattern (IBD; Orsini et al., 2013;Wright, 1943), where genetic distance increases with increasing geographic distance, but this scenario did not explain L. ocellata population structure as the correlation between genetic and geographic distances was not significant.
Vicariance, the process by which historical gene flow is disrupted by physical barriers such as mountain ranges, is one of the leading hypotheses to explain current distribution and phylogeographic structure of terrestrial animals (Brown et al., 2002;Pyron & Burbrink, 2010;Wiley, 1988). As L. ocellata were mostly confined to lowlands among mountain ranges in northern Thailand, we first considered that mountain ranges may play a significant role in their population structuring. However, this hypothesis was rejected due to the fact that the same haplogroups were observed on different sides of the mountains. For example, haplogroup C was present on both western (KY and MC), and eastern (CM, and CS) sides of the highest mountain of Thailand, Thanon Thong Chai range.
Furthermore, the Pleistocene divergence time was much later than the time of orogenic events in northern Thailand, which took place during the Tertiary (Tapponnier et al., 2001;Yin & Harrison, 2000).

This pattern was different from that of L. reevesii populations on
Hainan Island, where genetic differentiation occurred due to the presence of geographic barrier, the Wushishan mountain range (Lin et al., 2010), though elevation of the Wushishan range (1,840 m) was lower than that of the maximum elevation of the Thanon Thong Chai range (2,526 m).
River systems are much more likely to have act as geographic barriers to terrestrial species in the Indo-Burma region (Bain & Hurley, 2011;Geissler et al., 2015;Gonçalves et al., 2018;Klabacka et al., 2020;Rakotoarisoa et al., 2013). The rivers in northern Thailand, like the mountains, are oriented north-south.
As we have indicated, the Salween in the west and the Mekong in the east joined and flowed to the Siam river, presently recognized as the Chao Phraya river, during the Pleistocene (Hutchison, 1989;Woodruff, 2010). The capture of these river systems would have increased water volume such that they presented geographic barriers influencing the phylogeographic pattern of L. ocellata in northern Thailand (Klabacka et al., 2020). The Mekong river was thought to flow through the presently recognized Yom river (Brookfield, 1998) and may be responsible for the early divergence of L. ocellata. The ancestors of Leiolepis originated in the southern hemisphere and arrived in Southeast Asia approximately 120 mya (Macey et al., 2000), an ancestor of L. ocellata may have F I G U R E 4 Bayesian Inference (BI) tree based on Cyt b showed divergence time estimation of L. ocellata and L. belliana serving as an outgroup, using uncorrelated lognormal relaxed clock in BEAST v1.10.4. Filled circles with roman numbers indicated nodes of major divergences (details were in Table 3). Blue bars represented the 95% HPD (Bayesian credible interval) for timing of divergence F I G U R E 5 Predicted distribution of L. ocellata based on LGM, Holocene, and Present bioclimatic data. The highest probability of occurrence was shown in red color, while the lowest probability of occurrence was shown in blue color come from the south and distributed in southernmost lowlands in northern Thailand. The Yom river might have fragmented the ancestral L. ocellata into two major groups corresponding with node I in Figure 4. In eastern side of the lower Yom river, haplogroup A, principally located in DK, potentially diverged earlier, accounting for the relative high differentiation between the DK population and those of other localities. Part of the haplogroup A, population from DK, may have subsequently dispersed to adjacent localities such as WS, and also diverged as a closely related haplogroup B in MW (node II; Figure 4). The Salween river was hypothesized to flow through the presently recognized Ping river (Woodruff, 2010) in the west of northern Thailand. This may have isolated populations of L. ocellata found to the west of the Yom river into haplogroups C and D-F, mostly distributed in the west and central of northern Thailand, respectively (node III; Figure 4). This riverine vicariance due to the Ping river was supported by relatively high genetic differentiation between the populations on the west bank (KY and MC) and those populations to the east.
The phylogeographic pattern of L. ocellata was explained both by vicariance and also by dispersal following riverine divergence, as indicated by ecological niche modeling. The models revealed that distribution range of L. ocellata was largest during the LGM.
This was probably due to the fact that L. ocellata prefers dry habitat, which would have reached its maximum extent throughout northern Thailand at that time (Heaney, 1991;Penny, 2001 Such a dispersal scenario could explain the lineage diversification and admixture among haplogroups in the localities, located in the middle of the region between the Ping and Yom rivers (e.g., those at BT, MW, and TK). Moreover, dispersal of L. ocellata during the range expansions or minor changes in river course may have resulted in secondary contact between haplogroups confined to different banks of the rivers (Canestrelli et al., 2007;Ding et al., 2011;Graham et al., 2013). Localities CS and CM, for example, had an admixture of haplogroup C and F, which were found in the western and eastern banks of the Ping river, respectively. When the Salween river diverged from the Ping river, into its presently recognized course in Myanmar (Brookfield, 1998;Hutchison, 1989 (Brookfield, 1998;Hutchison, 1989). However, the presence of dominant haplogroup, instead of genetic admixture, in some localities (e.g., MN) may be a result of random effects such as genetic drift (Frankham, 1995).
The hypothesis of recent dispersal after riverine divergence is also supported by relatively small genetic differentiation among localities, except for KY, MC, and DK that diverged earlier (Kim et al., 2020).
In conclusion, the lineage diversification among L. ocellata populations has taken place too recently to have been influenced by the orogeny of the mountain ranges, but coincided temporally and spatially with the development of riverine geographic barriers.
These caused the early phylogenetic split of eastern and western populations. The original lineage diversification was then modified by subsequent dispersal and secondary contact followed by genetic admixture, especially in the central of northern Thailand. Our conclusions were derived from an integrative approach, combining population genetics, phylogeography, and ecological niche modeling, which enabled us to suggest hypotheses concerning the evolutionary history of L. ocellata in response to natural drivers in northern Thailand. We believed that further investigations of other species will assist clarification of phylogeographic puzzles in the northern Thailand biodiversity hotspot.

ACK N OWLED G EM ENTS
We are grateful to Jesse Grismer and his colleagues for providing some sequences. We would like to thank all those who assisted

CO N FLI C T O F I NTE R E S T
We have no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
All sequencing data are available at the GenBank, accession no.