The Silk roads: phylogeography of Central Asian dice snakes (Serpentes: Natricidae) shaped by rivers in deserts and mountain valleys

Abstract

Influenced by rapid changes in climate and landscape features since the Miocene, widely distributed species provide suitable models to study the environmental impact on their evolution and current genetic diversity. The dice snake Natrix tessellata, widely distributed in the Western Palearctic is one such species. We aimed to resolve a detailed phylogeography of N. tessellata with a focus on the Central Asian clade with 4 and the Anatolia clade with 3 mitochondrial lineages, trace their origin, and correlate the environmental changes that affected their distribution through time. The expected time of divergence of both clades began at 3.7 Mya in the Pliocene, reaching lineage differentiation approximately 1 million years later. The genetic diversity in both clades is rich, suggesting different ancestral areas, glacial refugia, demographic changes, and colonization routes. The Caspian lineage is the most widespread lineage in Central Asia, distributed around the Caspian Sea and reaching the foothills of the Hindu Kush Mountains in Afghanistan, and Eastern European lowlands in the west. Its distribution is limited by deserts, mountains, and cold steppe environments. Similarly, Kazakhstan and Uzbekistan lineages followed the Amu Darya and the Syr Darya water systems in Central Asia, with ranges delimited by the large Kyzylkum and Karakum deserts. On the western side, there are several lineages within the Anatolia clade that converged in the central part of the peninsula with 2 being endemic to Western Asia. The distribution of both main clades was affected by expansion from their Pleistocene glacial refugia around the Caspian Sea and in the valleys of Central Asia as well as by environmental changes, mostly through aridification.

Keywords: Eurasia; Paratethys; biogeography; colonization; genetic diversity; mitochondrial DNA; refugia; water snakes.

Figure 1

The maximum likelihood tree reconstruction of Natrix tessellata (1,117 bp; Supplementary Table S1) with time divergences based on the molecular clock analysis (A), and the inset to relationships between Anatolia and Central Asia clades (B). Numbers at nodes represent the expected time of the divergence. Terminal branch labels consist of the sample ID (new material) or GenBank accession number, and the name of the country or region of origin (see also Figure 1 and Supplementary Table S1). (C) Bayesian skyline plots showing the historical demography for the investigated lineages representing enough number of sequence data: the central line shows the mean value of the population size (N_e × τ × μ; where N_e is the effective population size, τ is the generation length in units of time [substitutions/site], and μ is the mutation rate) on the logarithmic scale. LGM line indicates the time of the Last Glacial Maximum. Inset photographs: Daniel Jablonski and M. M. Beskaravaynyi.

Figure 2

Upper panel (A): Geographic origin of sequences of the Anatolia and Central Asia clades of Natrix tessellata used in our study. Colors correspond to the main phylogenetic lineages recovered in our analysis (Figure 1; for locality details see Supplementary Table S1). Neighboring clades (Iran, Jordan, Europe, sensu Guicking et al. 2009) are indicated by different symbols. The distribution range of the species is highlighted in light orange. The question mark denotes uncertain origin of sequence KY887502. The pictured individual originates from Kazarman, Kyrgyzstan (Kazakhstan lineage). Lower panel (B): Ancestral areas of the genetic lineages of N. tessellata from Anatolia and Central Asia. Polygons represent regions with 10–70% highest posterior density (HPD) of the ancestral areas. The hypothetical colonization routes Caspian, Kazakhstan, and Uzbekistan lineage indicating by color arrows. Question marks indicate unknown genetic affiliation in unstudied areas and white arrows hypothetical origins. The map was drawn using QGIS 3.20. (https://qgis.org). Inset photograph: Daniel Jablonski.

Figure 3

The environmental niche model projection and the climate heterogeneity raster based on WorldClim data. (A) Climatic conditions at the Last Glacial Maximum (LGM). The white arrows suggest examples of potential glacial refugia of the species during the LGM. (B) The model for the Mid-Holocene. (C) Present model with the layer of data distribution points (white circles) downloaded from GBIF (2022). Yellow circles represent data used for the SDM projection. (D) Warm colors depict high areas of climatic heterogeneity in the present time. (E) Principal components analysis (PCA) of filtered WorldClim variables, showing climate space: the more similar the colors the more similar values. The maps were designed in QGIS 3.20 using Min–Max as the stretching histogram.