We used genetic, morphological, and environmental data to investigate diversity in rice paddy snakes (Homalopsidae, Hypsiscopus) and test if their diversification coincided with the Khorat Plateau uplift in Thailand. Previous studies on homalopsid snakes have found similar phylogenetic and morphological patterns[11, 27]. In this study, we increased genetic and geographic sampling density, and examined alternative population histories under a coalescent simulation framework. Multiple data types identified previously undocumented patterns of interspecific and intraspecific diversity in this genus. Our results suggest both geological processes and environmental factors have led to the present-day patterns of species diversity and their respective distributions.
Molecular and Morphological Diversity of Hypsiscopus
The concatenated phylogeny constructed with phased genetic data recovered each species as monophyletic. This topology is consistent with evolutionary relationships estimated previously, using both multilocus and genomic datasets[11, 27]. Our study includes newly added samples of H. matannensis from Sulawesi and H. murphyi near and east of the Red River in northern Vietnam (called the Red River population hereafter). The Red River population is monophyletic and sister to all other H. murphyi (Fig. 2). The Red River is considered to be a major biogeographic barrier[38], with known divergences occurring here for birds[82], frogs[38, 83–85], and lizards[86]. Within this Red River population, we observe substructure between samples from Hainan + Guangxi and those in Vietnam. One sample from the Vietnam clade is found on Hainan, which is not surprising given the land bridge between mainland China and Hainan Island during sea level minima in the Pleistocene and the close proximity of the Leizhou Peninsula to Hainan (~ 20 km).
Although the new samples of H. matannensis revealed substructure in the concatenated phylogeny (Fig. 2), the clades recovered did not reflect what may be expected based on Sulawesi’s geological history [87, 88]. Sulawesi’s arms (peninsulas) that branch from the island’s center (the Central Core) contain multiple potential contact points of paleo-islands when Sulawesi was partially submerged over the last 20 million years. Many of the inter-arm bays shrank in the last 4 million years leading to the coalescence of paleo-islands of Sulawesi, especially during the Pleistocene[88–91]. Many vertebrates show geographic-based population structure that reflects this history, such as Draco lizards[92], Bunomys rodents[93], Limnonectes frogs[94], and studies on macaque monkeys and Celebes toads have led to the establishment of areas of endemism (AOEs[95–98]), which are not reflected in Sulawesi rice paddy snakes. Some of our samples of H. matannensis from multiple peninsulas or islands showed minimal, intraspecific divergence and are recovered in multiple clades (Supplementary Fig. S10). Although further work and sampling is needed to thoroughly assess the biogeography of rice paddy snakes on Sulawesi, studies on lizards[99], shrews[100], and snails[101] have also found distributions that are partially incongruent with the monkey and toad AOEs or known historical geographic units (e.g., paleo-islands). This might be expected for more recent colonizations of Sulawesi that happened after the merger of paleo-islands had already taken place. Interestingly, the lasts paleo-island to merge was the southwest peninsula, which is the only place we detect any phylogenetic structure.
Hypsiscopus are absent from Australia, New Guinea, the Lesser Sunda Islands, and the Philippines. Thus, colonization of Sulawesi from the north, south, or east is unlikely. The Makassar Strait that separates Sulawesi and its west-neighboring island Borneo forms a strong barrier for reptiles and amphibians, with ~ 100 genera from Borneo being absent from Sulawesi[92]. Vicariance between Sulawesi and Borneo is ruled out, as this deep straight opened ~ 45 mya[102], significantly older than the crown date of Hypsiscopus. Thus, an over-water dispersal might be the most likely dispersal scenario into Sulawesi. Many species, such as terrestrial Sphenomorphus skinks[103], Draco lizards[104], Cyrtodactylus geckos[105],, and frogs in the genus Limonectes[94] have been suggested to have dispersed between islands on the Lesser Sunda Arc (Lombok, Sumbawa, Komodo, Flores, and Lembata). However, although our phased concatenated tree supports a Southwest colonization of Sulawesi (likely from Selayar Island, for which samples here cluster with those of the Southwest Peninsula), a southern dispersal from the Lesser Sunda Islands is considered improbable pending the discovery of Hypsiscopus from these islands. Future studies that include samples from the Central Core of Sulawesi and greater sampling of the Greater Sunda Islands will likely identify colonization routes into Sulawesi.
Phenotypic similarity is expected to be greater within species than between species. However, underlying species divergence is sometimes masked by conserved or parallel phenotypes in different taxa with wide distributions[23, 106, 107]. The morphological differentiation of H. plumbea and both populations of H. murphyi is quantitatively supported by the LDA (Fig. 2), with nearly perfect group classification accuracy (see Results: Morphological Quantitative Statistics). The Red River population of H. murphyi is phenotypically distinct from its conspecifics, having a gradual color change from dorsal to ventral scales and half-moon patterns on each ventral scale. This Red River population of H. murphyi is more phenotypically similar to H. plumbea than to other H. murphyi. It is unclear if the conserved color pattern in Red River H. murphyi and H. plumbea is due to incomplete lineage sorting of alleles related to color pattern or due to convergence in the environmental niches they occupy, as these species exist in both mainland (Indochina and Malay Peninsula) and island (Hainan, Taiwan, Greater Sundas) localities. Research on microhylid frogs also found in this region has shown ecological differences between the east and west sides of the Red River[38], though, our niche models show a difference in elevation suitability between H. murphyi (both populations) and H. plumbea, and the Red River population of H. murphyi is at an elevation more similar to the distribution of H. plumbea (see below). Both niche conservatism and niche divergence could drive the observable patterns seen in species. Niche conservatism can drive phenotypic conservatism, and both niche conservatism[108, 109] and niche divergence[108, 110, 111].
Geological and Environmental Drivers of Diversity
Our results support the hypothesis that Hypsiscopus diversified ~ 2.5 mya, suggesting the Khorat Plateau influenced the diversification of rice paddy snakes in Southeast Asia. All homalopsid genera most closely related to Hypsiscopus are restricted to mainland Southeast Asia[1, 11]. Demographic models showed no evidence of migration or change in Ne; this, and the distribution of outgroups, suggests that Hypsiscopus diversified from mainland Southeast Asia, eastward through the Greater Sunda Islands, eventually colonizing Sulawesi.
Pleistocene sea level fluctuations are one of the prevailing paradigms of diversification of land vertebrates—via sea level vicariance—in Southeast Asia[112, 113]. Phylogeographic patterns that temporally coincide with sea-level minima and maxima have been observed in spiders[114], freshwater fish[115, 116], birds[117], frogs[118, 119], gekkonid lizards[120], and shrews[121]. However, inland geological influences are a more likely explanation for some taxa. These include paleo-rivers and drainage basins[38, 86, 122, 123], mountains[124], and Quaternary landscape dynamics[125]. Tectonic uplift events, such as those that formed the Khorat Plateau, can also affect the distributions for elevation-sensitive species.
Although our BSP suggests changes in Ne may have occurred in the Pleistocene for each species (Supplementary Fig. S1), our demographic model results suggest that the much older Khorat Plateau uplift event might have initiated the divergence of rice paddy snakes in Southeast Asia and that no demographic changes have taken place (or these signatures were not strong enough to be detected). We note that our BSP were estimated using only mitochondrial data, but the demographic models take information from all loci (nuclear and mitochondrial), thus we favor the results of the demographic models that support no change in Ne over time. The Khorat Plateau is not particularly high in elevation (and geologically, is actually a basin with a depression in the middle, but we refer to it as a plateau due to its formal name in the literature). It consists of the Lower Mekong River, Phanom Dong Rak Mountains, and Phetchabun Mountains on the eastern, southern, and western sides, respectively, which reach elevations of 200–1,100 m above sea level. Even though the Khorat Plateau rim sites are not particularly high in elevation, homalopsids are generally considered a low-elevation group (Murphy, 2007). Additionally, studies on the Khorat Plateau show differing densities of some mud snake species (Enhydris spp., Subsessor bocourti, Erpeton tentaculatum, Homalopsis buccata, Hypsiscopus plumbea [sensu lato]) on the plateau in comparison to outside of the plateau, and no mud snakes of any species were found on the rim sites[126]. Furthermore, morphological differences in size and extent of sexual size dimorphism are significantly different when comparing populations on and off the plateau[126]. It is clear that the plateau leads to at least some inhibition of migration of mud snakes—a general result which has been seen in other vertebrates, including volant species like birds of the genus Alophoixus[127, 128], emydid turtles (Malalemys[129]), and potentially in Figs. 130]. It is worth noting that the basin is tilted towards the south and the east with an average elevation of 200 m, and the rim sites act as physical barriers that weaken the effects of the monsoons that commonly impact this region. This makes the plateau the hottest and driest area in Thailand with the greatest seasonal flux between wet and dry seasons and contains low-fertility soils[39, 131]. Thus, the Khorat Plateau not only represents a physical barrier, but may also serve as an ecological barrier between species as well. Environmental differences may affect the distributions of rice paddy snakes, and we see evidence of this in our niche modeling and quantification as well (Fig. 3).
Our ecological niche models reflect the currently known distributions of rice paddy snakes and suggest that their present-day distributions occupy different niches, most notably in elevation. H. murphyi showed higher suitability in moderate elevation areas, including the Khorat Plateau, with lower suitability at very high or very low elevations and south of the Isthmus of Kra biogeographic barrier on the Malay Peninsula. In contrast, niche models showed higher suitability at very low elevations for H. plumbea and H. matannensis (except on Sulawesi for H. matannensis), with the lowest suitability for both species in the chain of volcanoes that span the Greater and Lesser Sunda Islands and the high mountains of central Borneo (~ 1,000–4,000 m elevation). This suggests environmental differences in elevation, as well as annual mean temperature, mean diurnal range, temperature seasonality, mean temperature of warmest quarter, annual precipitation, precipitation seasonality, precipitation of wettest quarter, and precipitation of warmest quarter (Supplementary Fig. S3) likely impact the current distribution of Hypsiscopus. The ENM projections to the LGM showed differences in habitat suitability between all species, suggesting that environmental differences may be responsible for historical diversification processes leading to the extant species richness of these snakes. Our niche quantification analysis supports the interpretation of niche differentiation between H. murphyi and H. matannensis, reinforcing the ecological differentiation observed between these closely related taxa, possibly via a gradual evolutionary niche shift in Hysiscopus as the ancestors of today’s lineages diversified from west to east (Fig. 4).
Our work suggests that the diversification of Hypsiscopus may have been influenced by multiple processes, with an initial vicariant divergence due to the Khorat Plateau uplift, and then distinct environmental features (both near and far from the plateau) helped maintain their divergence. Environmental differences in mainland Southeast Asia have led to the diversification of microhylid frogs[38] and birds[132, 133], the latter of which many have near-identical distributions as H. plumbea. We highlight the importance of utilizing multilocus datasets and numerous methods and data types to test fine-scale biogeographic hypotheses and identify patterns of differentiation in poorly studied, yet widespread species. The evolutionary history of rice paddy snakes has only recently been studied[27]. Our sampling of nearly 200 individuals provides a basis for future studies to investigate broad and fine-scale processes such as gene flow, and dispersal—and to characterize the colonization routes both may take throughout the complex geography of this environmentally heterogeneous region. Our empirical approach in this study allows us to ultimately understand dynamic geological and environmental processes that continue to shape Southeast Asia, as well as the biodiversity generated, partitioned, and maintained on either side of Wallace’s Line.