The biogeography of bent-toed geckos, Cyrtodactylus (Squamata: Gekkonidae)

Biogeographic regions (Fig. 1)

A presence or absence species geography file was constructed using generally the same biogeographic areas employed by Wood Jr et al. (2012) with modifications and further refinement herein (see below). These areas were not randomly chosen but are regions among which cladogenic turnover within a broad range of taxonomic groups has been demonstrated. No species ranged outside their designated biogeographic region and as such was allowed to occupy only a single area.

Overview of regional dispersal patterns of Cyrtodactylus (Fig. 4)

Cyrtodactylus evolved during the early Eocene approximately 52 million years ago (mya) and gave rise to two small monophyletic Himalayan lineages before dispersing into what is currently Indoburma and Indochina. This is consistent with the step-wise scenario proposed for the Himalayan region by Xu et al. (2020). The colonization of Indoburma resulted in the evolution of a modest number of species but that of Indochina seeded the origin and evolution of the remainder of the genus from India and Sri Lanka to Melanesia (see alsoWood Jr et al., 2012 and Xu et al. (2020). Three independent, consecutive, episodes of dispersal out of Indochina gave rise to radiations in the regions of Sundaland, India-Sri Lanka, and Papua. The Sundaland radiation spread eastward and gave rise to groups in Wallacea on at least three separate occasions with dispersals back and forth between the two regions for approximately 15 million years, during which there was also a dispersal to Papua from Wallacea. The Sundaland radiation also gave rise to a Philippine clade and a single dispersal back to Indochina with a subsequent radiation therein. The India-Sri Lanka radiation spread across southern and central India and onto Sri Lanka. The Paupa radiation gave rise to several small, but distinctive clades throughout that region, as well as one clade that dispersed back west to colonize and radiate in Sundaland. Currently, distantly related and independently evolved clades throughout Indochina and Sundaland account for the majority of the diversity within Cyrtodactylus and continue to grow at an ever-increasing rate with ongoing discoveries of new species (see Uetz, Freed & Hošek, 2021).

Evolution, dispersal, and colonization of Cyrtodactylus (Figs. 1 and 3)

Below we discuss the biogeographic patterns of the major clades and species groups recently delimited in Cyrtodactylus (Grismer et al., 2021a; Grismer et al., 2021b). We do not discuss many of the fine-scaled biogeographic patterns within these clades and species groups that are detailed elsewhere, but instead, direct the reader to these fine publications at the relevant junctures. All ages presented are mean values and the highest posterior densities (HPD) about the means are presented in Fig. 2. Wood Jr et al. (2012) and Agarwal et al. (2014) noted that a number of cladogenic events among basal lineages of Cyrtodactylus and other gekkonid genera occurred early on the southern edge of the Asian continent in a “proto-Himalayan” region, or near what is currently the Tibetan Plateau from ∼62–52 million years ago (mya). These divergences were hypothesized to result from changing landscapes during ongoing compression of the Neotethyan island arcs, including the Kohistan Arc—situated between the Indian subcontinent and Eurasia (Ali & Aitchison, 2008; Morley, 2018)—as the Indian subcontinent began approaching continental Asia (Fig. 5A). These events are proposed to have sequentially given rise to the ancestor of the tibetanus group at ∼51.6 mya, and the ancestor of the lawderanus group, and ancestor (A)1 ∼47.7 mya (Fig. 3). The dates here are slightly younger (∼5 million years) than those of Agarwal et al. (2014) but are commensurate with those of Wood Jr et al. (2012) upon which our tree was calibrated.

Ancestor A1 diverged in the Himalayan region ∼43.8 mya, giving rise to A2—which remained in the Himalayan region—and A3. A3 dispersed to the Indochina region through mountainous corridors along the southeastern edge of an emerging Qinghai-Tibetan Plateau (Farve et al., 2015) as these basal divergences formed chronologically separate lineages in a west to east pattern from the Himalayas to Indoburma (Wood Jr et al., 2012; Agarwal et al., 2014). One of these divergences involving A2 in the Himalaya region, occurred ∼39.3 mya and resulted in the formation of the fasciolatus group and A4—the latter of which dispersed into the lowlands of the Indoburma region along north-south oriented hilly corridors (Farve et al., 2015; Ding et al., 2017) where it radiated ∼35.4 mya into the peguensis and khasiensis groups. This coincides with the intensifying orogeny in the region, which during the early Oligocene formed the main Indoburman ranges: Arakan, Chin, and Naga Hills (Brunnschweiler, 1966) and would have facilitated the invasion of new habitats and speciation from a habitat generalist ancestor (Grismer et al., 2020). The peguensis group began to radiate ∼33.2 mya from the southern foothills of the Himalayas to the eastern edge of the Ayeyarwady Basin (Grismer et al., 2021a: Fig. 24) and the khasiensis group began to radiate ∼29.5 mya from generally south of the Brahmaputra River and east of the Padma River to the eastern edge of Ayeyarwady Basin (Agarwal et al., 2014; Grismer et al., 2021a: Fig. 15). The biogeography of the species within the khasiensis group is discussed further in Agarwal et al. (2014). The fasciolatus, khasiensis, and peguensis groups were referred to as the Myanmar clade in Wood Jr et al. (2012).

Orogenic events in what is currently northern Indochina between ∼40 and 35 mya, may have precipitated sequential, cladogeneic events among the Indochinese ancestors A3, A6, A12, and A14. A3 diverged into two lineages within Indochina ∼40.9 mya. One lineage, A5, radiated in karstic habitats throughout central Indochina from approximately 33.1 mya to present, giving rise to the angularis group (Grismer et al., 2021a: Fig. 9). The other lineage A6, diverged ∼39.5 mya and gave rise to A12, which remained in Indochina, and A7 which dispersed into Sundaland—most likely into current-day northern Borneo as they were a contiguous landmass during this period (Hall, 2013; De Bruyn et al., 2014; Fig. 5; Cao et al., 2017). This portion of Sundaland remained subaerial throughout the tectonic evolution of Southeast Asia, was covered with perhumid rainforests since at least the middle Eocene, and became a major center of origin for a vast number of clades of plants and animals (De Bruyn et al., 2014; Grismer et al., 2016; Morley, 2018, and references therein), including Cyrtodactylus (see below).

As A7 began to diverge ∼33.1 mya in Sundaland it became the ultimate source of origin for many species and species groups throughout Sundaland, Wallacea, Papua, and the Philippines, and included one back-dispersal into Indochina. The cladogenic events within the descendants A7 were likely the result of the dynamic and rapidly fluctuating insular landscapes of Sundaland from ∼30–10 mya (Hall, 2013; Fig. 5). A7 gave rise to a clade comprised of the lateralis and sworderi sister groups which radiated in western Sundaland (Grismer et al., 2021a: Fig. 20) from ∼27.7 mya to present. Grismer & Davis (2018) and O’Connell et al. (2019) discuss the biogeography of the laterlais and sworderi groups. A7 eventually gave rise to three ancestors that independently dispersed into Wallacea. One, presumably made an over-water dispersal onto the western Sulawesi landmass sometime after ∼28.2 mya when it diverged from the malayanus group and is represented today by C. spinosus. The western part of Sulawesi was connected to the Sunda shelf until at least 25 mya (Cao et al., 2017). The malayanus group radiated in Borneo from ∼20.2 mya to present. Davis et al. (2020) discuss the biogeography of the malayanus group. Wood Jr et al. (2012) hypothesized that regional endemism in mainland Southeast Asia may have been facilitated by long-term geographic barriers and is further supported here with substantially more species sampling (see below).

Another descendent of A7, A8, dispersed into Wallacea ∼29.1 mya, giving rise to the darmandvillei group which began to radiate ∼23.4 mya. Members of the darmandvillei group dispersed over-water back and forth across what is currently recognized as Wallace’s Line into Sundaland (O’Connell et al., 2019) at least three times between ∼19.6–5.3 mya (Grismer et al., 2021a: Fig. 7). The colonization of northwestern Australia by C. kimberleyensis happened no earlier than 5.3 mya when it diverged from its undescribed sister species on Timor Island (Wood Jr et al., 2012). Cyrtodactylus sadleiri on Christmas Island diverged no earlier than 3.9 mya from its undescribed sister species on Bali Island, Indonesia. A third Sundaic descendent, A9, diverged ∼29.1 mya and gave rise to a clade comprised of the agamensis and marmoratus groups, likely as a result of intensifying volcanism and orogeny on the western edge of Sundaland (Hall, 2012). These groups separated in Sundaland ∼27.7 mya and the agamensis group radiated in western Sundaland from approximately 24.3 mya to present (Grismer et al., 2021a: Fig. 7), whereas the marmoratus group began radiating in southern Sundaland ∼23.7 mya and dispersed to the Papua region through Wallacea (Grismer et al., 2021a: Fig. 27) ∼10.4 mya. Grismer & Davis (2018) and O’Connell et al. (2019) discuss the biogeography of the agamensis group.

As sequential cladogenic events—first initiated in A7—continued throughout Sundaland, most likely in what is current-day Borneo, one lineage ∼29.9 mya, gave rise to ancestor A10 of the philippinicus group which sequentially gave rise to two Bornean clades (Wood Jr et al., 2012). The first at ∼26.7 mya, is currently represented in part by the Cyrtodactylus pubisulcus species complex (Davis et al., 2020). The second Bornean clade, represented today by at least 10 species, diverged ∼24.3 mya and radiated in Borneo while its sister lineage dispersed to the Philippines sometime afterwards. The biogeography of the species of both Bornean clades are discussed by Davis et al. (2020). The Philippine clade began radiating in situ ∼23.1 mya, indicating its route to this archipelago was along a southern chain of islands onto Mindanao as it was the only available route during this time period (Hall, 2013; Fig. 5). Siler et al. (2010), Siler et al. (2012), Welton et al. (2010a) and Welton et al. (2010b) discuss the biogeography of the Philippine gekkonids and other possible routes. Eventually, A7 gave rise to an Indochinese descendent A11 ∼30.1 mya, which dispersed back into the Indochina region no later than ∼27.5 mya during a time period when broad, subaerial upland areas still linked northern Borneo and southern Indochina across an exposed Sunda Shelf (Figs. 5D–5E). A11 diverged, giving rise to a clade of at least 33 species comprising the condorensis and irregularis groups. The condorensis group is centered in southern Vietnam and radiated from ∼23.7 mya to present and the irregularis group radiated in Vietnam, Laos, and Cambodia from ∼23.3 mya to present (Grismer et al., 2021a: Fig. 9). One species of the condorensis group, C. leegrismeri, occurs on two small islands off the east coast of the Burma-Thai-Malay Peninsula in Sundaland. Grismer & Grismer (2017) discuss the biogeography of the condorensis group and Nguyen et al. (2017) discuss the biogeography of the irregularis group.

Approximately 38.1 mya, ancestor A12 diverged in Indochina, giving rise to the Indochinese ancestor A14 and the ancestor of the triedrus group (A13) which dispersed to the India-Sri Lanka region at a time when the Indian subcontinent was adjacent to Indoburma and Indochina (Acton, 1999; Köhler & Glaubrecht, 2007; Aitchison, Ali & Davis, 2007; Ali & Aitchison, 2008: Fig. 5). This invasion route to Indian-Sri Lanka from what was probably southern Indochina (being there is no evidence of this group ever occupying Indoburma), via a relatively narrow over-water dispersal or land bridge, has been hypothesized for a number taxonomic groups (e.g., (Klaus et al., 2010) [crabs]; (Li et al., 2013) [rhacophorid frogs]; J. (Grismer et al., 2016) [draconine lizards]; (Garg & Biju, 2019) [microhylid frogs]; (Gorin et al., 2020) [microhylid frogs]). Given the extensively greater sampling of Cyrtodactylus here and the increased use of tectonic data, this scenario does not support long distance over-water dispersal scenarios across the Bay of Bengal proposed by Wood Jr et al. (2012) and Agarwal et al. (2014) nor does it support the scenario of Grismer & Davis (2018) that Sundaland was colonized by an ancestor of the Indian subcontinent. Based on the above interpretations, the triedrus group began diversifying in mesic areas no later than ∼33.3 mya and has become generally restricted to the periphery of the subcontinent south of the Ganga River following gradual aridification during the Oligocene (Morley, 2018; Deepak & Karanth, 2018). Agarwal & Karanth (2015) discuss the biogeography of the terrestrial members of this group.

A14 remained in Indochina and continued to radiate, eventually diverging ∼35.8 mya and giving rise to A19, which also remained in Indochina, and A15 which invaded the Papua region. This colonization event is problematic. Unlike the invasion of the Papua region by the marmortus group which left a chronological sequential trail of ancestors throughout Sundaland and Wallacea, no such trail exists for A15. From ∼40–30 mya, a series of islands that later coalesced to form Sulawesi existed to the east of southern Indochina—which at the time was contiguous with current-day Borneo—and could have provided a route into the Papua region (Hall, 2013; Fig. 5). The absence of closely related species in current-day Sundaland and Wallacea neither supports nor discounts this hypothesis. A15 diverged ∼32.7 mya in the Papua region and gave rise to A16—the ancestor of a monophyletic Melanesia radiation that currently consists of at least 36 species comprising seven monophyletic species groups that radiated throughout the Papua region (including northeastern Queensland) from ∼29.5 mya to present (Grismer et al., 2021a: Figs. 37, 39). Oliver, Skipwith & Lee (2014) and Oliver et al. (2018) discuss the complex biogeographic history of this clade. At 32.7 mya, A15 also gave rise to A17 in Papua New Guinea which diverged ∼28.2 mya and may have given rise to Cyrtodactylus biordinis in the Solomon Islands (Oliver et al., 2018). The putative sister lineage of C. biordinis, A18, dispersed to and radiated on the Burma-Thai-Malay Peninsula in western Sundaland south of the Isthmus of Kra as the pulchellus group (Grismer et al., 2021a: Fig. 16) from ∼18.3 mya to present. Like the dispersal from the Indochina region to the Papua region, this too is a difficult scenario to explain given that there are no closely related, earlier diverging species in the intervening regions of current-day Wallacea or in southern Sundaland. Grismer et al. (2014) discuss the biogeography of the pulchellus group.

A19, a direct descendent of A14 in the Indochina region, is the ancestor of a large Indochinese clade containing at least 78 species that comprises eight species groups. Grismer & Davis (2018) discuss the complex tectonic history that this vast, corrugated landscape of intermontane rift basins and parallel mountain ranges (Feton et al., 1997; Fenton, Charusiri & Wood, 2003) played in the diversification of this clade from eastern Myanmar to central Thailand and others (Chen et al., 2017; Chen et al., 2018; Grismer et al., 2018; Poyarkov et al., 2018; Suwannapoom et al., 2018; Gorin et al., 2020), as the Tenasserim Mountains and Thai Highlands were gradually uplifted (Upton, Bristow & Hurford, 1995; Upton et al., 1997). As A19 and its descendants continued to diverge and radiate in Indochina, the first ancestors to evolve were A20 and A21, ∼34.2 mya. A21 diverged approximately 29.7 mya, giving rise to the brevipalmatus group centered in Indochina and the linnwayensis group of the Shan Plateau of eastern Myanmar (Grismer et al., 2021a: Figs. 9, 24). Approximately 21.0 mya, the breviplamatus group dispersed across, what is currently the Isthmus of Kra on the Burma-Thai-Malay Peninsula, into northern Sundaland to as far south as Peninsular Malaysia where it is currently represented by C. elok and C. breviplamatus.

Ancestor A20 continued to radiate in the Indochina region and at ∼32.9 mya, it diverged, giving rise to A24 of the intermedius group which then radiated throughout eastern Thailand, Cambodia, and southern Vietnam from ∼24.6 mya to present (Grismer et al., 2021a: Fig. 16). A20 also gave rise to ancestor A22—a clade of Indochinese species comprised of the sadansinensis, yathepyanensis, oldhami, sinyineensis, and chauquangensis groups. The earliest diverging lineages of that clade are the sadansinensis (∼30.1 mya) and the yathepyanensis (∼28.0 mya) groups from the Salween Basin of southeastern Myanmar (Grismer et al., 2021a: Fig. 24). The next group to diverge, the oldhami group (∼26.7 mya), is endemic to the Burma-Thai-Malay Peninsula and contains one species, C. zebraicus, that diverged in northern Sundaland ∼13.4 mya, just south of the Isthmus of Kra (Grismer et al., 2021a: Fig. 16) and another undescribed Sundaic species from Krabi, Thailand, that is closely related to C. sanook, which also diverged south of the isthmus ∼5.5 mya. The most recently diverged sister lineages, the sinyineensis and chauquangensis groups, radiated throughout Indochina in allopatry with respect to one another. The sinyineensis group is a disjunctly distributed, saxicolous lineage that began radiating ∼22.3 mya and extends from the Salween Basin of southeastern Myanmar through the mountainous terrain of western Thailand to at least Doi Inthanon (Grismer et al., 2021a: Fig. 24). The chauquangensis group radiated from ∼21.7 mya to present, across karstic landscapes in northern Vietnam, southern China, Laos, and northern Thailand (Grismer et al., 2021a: Fig. 9).

Given that the mtDNA phylogeny here is highly congruent with the mito-nuclear phylogeny of Grismer et al. (2021a), we suspect that adding more species and genes will have little bearing on the phylogeny and the biographical scenario presented here. However, we are in the process of constructing a genomic data set with additional species that will be compared to the mitochondrial and mito-nuclear data sets of the most recent genus-wide phylogenies of Grismer et al. (2021a) and Grismer et al. (2021b).