What Makes a Hot-Spring Habitat “Hot” for the Hot-Spring Snake: Distributional Data and Niche Modelling for the Genus Thermophis (Serpentes, Colubridae)

: Knowledge about species’ distributions is central to diverse applications in ecology, biogeography, and conservation science. Hot-spring snakes of the genus Thermophis share a distribution restricted to geothermal sites at the Tibetan Plateau ( T. baileyi ) and in the Hengduan Mountains ( T. zhaoermii , T. shangrila ). Although the suture zones of these regions are widely covered with hot springs, Thermophis populations are restricted to only a few of these habitats. Here, we use bioclimatic, topographic, and land cover data to model the potential distribution of the genus. Moreover, using logistic regression on ﬁeld survey data of T. zhaoermii , we test whether hot-spring water parameters and landscape features correlate with the species’ presence or absence. Hot springs with temperatures between 45 and 100 ◦ C and winter precipitation showed the most predictive power. At small scale, our data support the relevance of the hot-spring temperature on the species’ occurrence and indicate that also the along-valley distance from the hot-spring site to the major river might inﬂuence the distribution of Thermophis species. Our ﬁndings contribute to better understand factors shaping the current distribution of the genus and will aid in setting priorities in applied conservation biology for the hot-spring snakes.


Introduction
Understanding the factors underlying the distribution patterns of species is one of the major tasks of ecology and biogeography [1,2]. It is, moreover, of prime concern for biodiversity conservation.
The Tibetan Plateau and the Hengduan Mountains pose as unique biogeographic zones since they are among the world's ecologically most diverse areas [3,4]. The Hengduan Mountains comprise more than 500,000 km 2 of temperate and alpine ecosystems that run along a north-south orientation at the eastern margin of the Himalayan range, in the southeastern corner of the Tibetan Plateau. Its southeastern quarter is one of the world's most important biodiversity hotspots characterized by high levels of species richness and endemism [5]. Similarly, the Tibetan Plateau supports high species diversity. With an area of 2.5 million km 2 and an average elevation of 4500 m above sea level (a.s.l.), the highlands Records are according to [17,18,21] and based on our survey data. Locations where we did not find Thermophis are included in the map (absence). Hot-spring sites are drawn according to the "Atlas of Tibet" [22]. Type localities of the three species are imprecise in the original descriptions and cannot be mapped exactly-T. baileyi: "Thibet"; T. shangrila: "Shangri-La, Northern Yunnan, China" but coordinates match Yulong, Lijiang, Yunnan Province, China; T. zaoermii: Litang County, Sichuan, China (probably locality no. 12 according to [20]). The extension of the Hengduan Mountains region (HMR) is roughly reflected by its typical temperate coniferous forests ecoregion [23]. Records are according to [17,18,21] and based on our survey data. Locations where we did not find Thermophis are included in the map (absence). Hot-spring sites are drawn according to the "Atlas of Tibet" [22]. Type localities of the three species are imprecise in the original descriptions and cannot be mapped exactly-T. baileyi: "Thibet"; T. shangrila: "Shangri-La, Northern Yunnan, China" but coordinates match Yulong, Lijiang, Yunnan Province, China; T. zaoermii: Litang County, Sichuan, China (probably locality no. 12 according to [20]). The extension of the Hengduan Mountains region (HMR) is roughly reflected by its typical temperate coniferous forests ecoregion [23].
Compared with the Tibetan Plateau where hot springs are much more sparsely distributed, the Hengduan Mountains are densely covered with hot springs [22] (Figure 1). Surprisingly, the distribution of T. zhaoermii and T. shangrila is virtually unknown for the Hengduan Mountains and adjacent mountain ranges, despite multiple herpetological field surveys that have been carried out in these regions (e.g., [24][25][26][27][28]). We suppose that this apparent contradiction between the large number of potential habitats and lack of hot-spring snake records might be explained by specific environmental characteristics that impact the suitability of a hot-spring locality for species of Thermophis and, thus, determine the presence or absence of this taxon in apparently similar habitats. Therefore, in the present study we aimed to model the habitat suitability and potential distribution of the genus to determine the environmental factors related to its occurrence. We further sought to identify chemical, physical, and spatial hot-spring water characteristics that may impact the distribution of the hot-spring snakes by collecting a comprehensive data set on water quality of surveyed hot springs. Our data represent the first evaluation of environmental features and occurrences of this enigmatic genus. The results encourage discussion on our understanding of the distribution of Thermophis and can have important implications for the management and conservation of these rare high-elevation snakes.

Occurrence Data and Ecological Niche Analyses
A total of 20 records of T. baileyi, 6 records of T. zhaoermii, and 1 record of T. shangrila were available to us (Figure 1, Supplementary Materials Figure S1). All of them, except the locality of T. shangrila, were obtained from field surveys during 2003-2014 [17,29] (Supplementary Text S1 and Table S1). Since many hot springs are located in remote, inaccessible areas, a systematic survey of such sites is difficult, if not impossible due to logistic and financial challenges. Therefore, locations were selected based on hot-spring literature references [22] and information provided by residents. Sites were visited between May and August, and were surveyed intensively for up to four days, depending on weather conditions (for details on site visits see Supplementary Text S1). We are aware that the absence of secretive species like snakes from a given site can only be determined with a certain probability and strongly depends, among other factors, on population density and number of visits [30]. However, hot-spring snakes are not difficult to find if looked for under suitable weather conditions during their active season. The availability of snake distribution data and the assistance of local residents were essential components of our successful searches. Geographical coordinates and altitude (ALT) were recorded with a GPS receiver using the WGS84 datum and imported into ArcGIS 10.8 (ESRI, Redlands, CA, USA). Overall, we surveyed 66 habitats for presence or absence of T. baileyi (n = 44; [17]) and T. zhaoermii (n = 22; Supplementary Table S1). The third species, T. shangrila, is known only from the holotype location and a further site in Shangri-La (Dêqên). However, since the coordinates given in the original description do not match the correspondingly specified locality "Shangri-La" [20] but Yulong, Lijiang, Yunnan Province, China, we excluded this record due to its spatial discrepancy. Due to the low number of reported localities for the three species and the general strong preference of the genus for hot-spring habitats, we combined all Thermophis records into a single data set.
Grids of 19 standard bioclimatic variables for the current climate and elevation were downloaded from the WorldClim database [31] (http://www.worldclim.org (accessed on 7 June 2021)). Hot springs were digitized from, and according to "The Atlas of Tibet" [22] using ArcGIS 10.8, subdividing them into the following groups (layers): 5-35 • C, >35-45 • C, >45-100 • C, and >100 • C. All layers were projected to WGS84 at a resolution of 30 arc-seconds (~1 km grid cells at the equator) and clipped to the extent of Tibet s.l. (incl. the northern part of the Indian state Arunachal Pradesh), and Sichuan and Yunnan Provinces, which cover all known Thermophis records, and which together correspond roughly to the main area where hot springs exist ( Figure 1). Prior to final model construction, we explored all climate variables and the elevation data for autocorrelation by calculating Pearson's correlation coefficients (r) using the python script SDMtoolbox v.2.4 [32] available for ArcGIS 10.8 and removed highly correlated variables (r > 0.8). The final set of environmental predictors comprised six variables: BIO1 = annual mean temperature, BIO2 = mean diurnal temperature range, BIO3 = isothermality, BIO4 = temperature seasonality, BIO15 = precipitation seasonality, and BIO19 = precipitation of coldest quarter. Principal components analysis (PCA) was carried out for these variables with the SDMtoolbox, reducing them to three orthogonal principal components describing the majority (>99%) of the variability in climate. We then used these three components to assess the climatic heterogeneity across the area of interest. To eliminate spatial clusters of localities we spatially filtered our presence data by Euclidian distances (min 2 km, max 20 km; 5 distance classes) according to climate heterogeneity using the rarefying module in SDMtoolbox. This resulted in the exclusion of one site (no. 15, see Supplementary Figure S1, Table S1). We also included land cover (LC) information in the model using GlobCover Land Cover V2.2 provided by the European Spatial Agency (http://due.esrin.esa.int/page_globcover.php (accessed on 7 June 2021)), as well as the river network (RN) using vmap0 data obtained from the National Geospatial-Intelligence Agency (https://www.nga.mil (accessed on 7 June 2021)). The original land cover map was rescaled to match the 30 arc-seconds resolution of the other variables. Vector data were converted to a raster using the same cell size as all other layers.
The final models were generated with the Maxent algorithm [33,34] which is extensively used for analyzing presence-only data and based on the principle of maximum entropy. Maxent is known for its high predictive accuracy compared with other modelling methods [35], and robustness to small sample sizes [36]. In the occurrence dataset, we implemented the filtered 25 georeferenced localities corresponding to the surveyed populations of Thermophis ( Figure 1, Supplementary Table S1; [17]). These occurrences were randomly split into training data (75%) and test data (25%), and 100 subsampled (model 1) or bootstrapped (model 2) replicates were run for model evaluation using the following parameters: regularization = 1.5; max. iterations = 50,000; min. threshold = 0.00001; output = logistic. Both of the two models comprised BIO1-4, 15, 19, LC, RN, and the hot-spring layers. Model performance and the importance of the environmental variables to the model were assessed using the mean area under the curve (AUC) of the receiver operating characteristics (ROC; [37,38]), and jack-knife testing. Jack-knife testing obtains models by first omitting each variable and then using only one variable to determine the importance of an environmental variable to the predictive distribution of a species [33]. The AUC ranges from 0 to 1 and can be interpreted as the probability that a randomly chosen presence site will be ranked above a randomly chosen absence or pseudo-absence (background) site [39]. An AUC value close to 1 indicates higher model fit, a value of 0.5 implies random prediction, and values <0.5 indicate performance worse than random; models with values above 0.75 are considered potentially informative [35], good between 0.8 and 0.9, and excellent for AUC between 0.9 and 1 [40]. Although, shortcomings of AUC have been made evident [41,42], we considered that criterion a reliable measurement to compare models generated for the same genus, within the same area and under the same settings.

Hot-Spring Water Characteristics
To include microhabitat variables in our assessment, in particular water characteristics of the hot springs, we additionally obtained several water chemistry measurements from a total of 16 hot-spring sites that were surveyed for T. zhaoermii in the Hengduan Mountains. All the surveyed sites, except for no. 1, were located in Sichuan Province, China (Supplementary Figure S1; Table S1) and presented the following characteristics, which were common to all habitats where the species has been recorded [17]: sites were in river valleys with rocky slopes and vegetated shorelines, within a short distance to a river (0-720 m), and could be considered as low-sulfur sites (i.e., without characteristic sulfur smell).
Since hot-spring snakes are suggested to spread mainly along river valleys [17], we assume increasing distance of a hot-spring habitat from a (main) riverine corridor to cause higher costs for the migratory pathway and, thus, to lower the probability that a Thermophis population has been established itself in that habitat. Therefore, we determined the length of the riverine corridor from the respective hot spring (located in the river valley) to the point where this river merges with the superordinate river and, further, with a major river that runs through the Hengduan Mountains (Dajin/Dadu, Jinsha, Litang or Yalong River) to explore the impact of the corridor length on species presence or absence. Data on distances of a hot-spring site to (i) the nearest small (seasonal) river (SRD), (ii) the next merging point with the superordinate river (tributary flow path of SRD; TRIB), and (iii) to the junction with the major river (MRD) were extracted as pathways along the river valleys based on Google Earth layer using ArcGIS 10.8. Measurements of hot-spring water chemistry were performed under field conditions next to the respective hot spring. We recorded the following parameters with a portable digital data logger PCE-PHD 1 (PCE Instruments, Meschede, Germany): temperature (TEMP; 0.1 ± 0.8 • C), conductivity (COND; 0.1 µS ± 2%), total dissolved solids (TDS; 1 ppm ± 2%), pH (0.01 ± 0.02 pH), oxidation reduction potential (ORP; 1 mV ± 0.5%), salinity (SAL; 0.01 ± 0.05%), and dissolved oxygen (OX; 0.1 ± 0.4 µg/L). Measurements of COND, pH, and OX were automatically temperature compensated. OX was moreover manually compensated for altitude. In order to determine several ionic concentrations, we used semi-quantitative (colorimetric) test kits (NO 3 − , S 2− , PO 4 3− , NH 3 /NH 4 , carbonate hardness: JBL GmbH, Neuhofen, Germany; S 2− : HC393108 and HC41474, Merck, KGaA, Darmstadt, Germany; SO 4 2 : Hanna Instruments, Vöhringen, Germany). Since the Maxent models showed the importance of BIO19 for the distribution of Thermophis, we extracted these data for the surveyed sites from WorldClim.

Logistic Regression Models
We used binomial regression models to test whether certain hot-spring water characteristics, climatic, and/or landscape features were favorable for the species' presence or absence (modelled as dependent variable and coded as 1 and 0, respectively). Given the very small sample size of our dataset, we were left to avoid model over-fitting of the data by applying exclusively an exploratory approach. In a first step, we calculated Spearman correlation coefficients for all pairs of water characteristics and removed variables with coefficients >0.7 in order to avoid multicollinearity (Supplementary Table S2). Next, we performed a forward model selection based on Akaike's information criterion (AIC; [43]), adding the remaining variables stepwise to the null model, to explore which variables could best discriminate between presence and absence of the species. All statistical analyses were carried out using the glm, step, and anova functions in R v.4.1.0 [44].

Climate Heterogeneity, Maxent Models
For the set of six analyzed climate variables, 99.7% of the variance was explained by the first ordination axis (Supplementary Figure S2A). Overall, climate across the modelled area is spatially heterogeneous. Areas in the Yarlung River catchment at the Tibetan Plateau and in the central Hengduan Mountains, where Thermophis occur, share similar climatic conditions (indicated by similar, blue-reddish color, see Supplementary Figure S2A), while regions between and around these areas are different. Particularly southeastern Tibet, with a topography dropping below 3000 m a.s.l. in many valleys, shows varying and higher mean annual temperatures and precipitation rates. Similarly, northern areas of Tibet and Sichuan, central and southern Yunnan, and the Sichuan Basin are climatically different from those areas where Thermophis has been recorded. Regions of higher climate heterogeneity can be found in the southeastern part of Tibet, the southern Hengduan Mountains, and along the margin of the Sichuan basin; the known distribution of Thermophis matches climatically more homogeneous areas (Supplementary Figure S2B).
The two Maxent models produced consistent results and had high predictive power (AUC m1 0.98 and AUC m2 0.99; Figure 2, Supplementary Figure S3). The test omission rate of m1 (subsample approach) is close to the predicted omission, while the omission rate Diversity 2021, 13, 325 6 of 11 of m2 (bootstrap approach) slightly deviates from the predictive values for cumulative thresholds <0.5 (Figure 2, Supplementary Figure S3).
Hot springs with a temperature between 45 • C and 100 • C and the amount of precipitation in the coldest quarter of the year (i.e., winter snowfall; BIO19) had the highest percentage contribution to the models (58%, 15% in m1, and 56%, 14% in m2, respectively; Figure 2, Supplementary Figure S3, Table S2). Both maps of the two models show very similar patchy distribution patterns of habitat suitability with high values (yellow to red) corresponding to pixels that overlap with several hot-spring sites in the central Hengduan Mountains and in the tributary valleys of the Yarlung River in southern central Tibet (Figure 2, Supplementary Figure S3). These suitable habitats are located in areas with lower winter precipitation than surrounding regions (Figure 3; Supplementary Figure S4). Hot springs outside these ranges are associated with low values of habitat suitability (green), and areas without hot springs are unsuitable (blue). heterogeneity can be found in the southeastern part of Tibet, the southern Hengduan Mountains, and along the margin of the Sichuan basin; the known distribution of Thermophis matches climatically more homogeneous areas (Supplementary Figure S2B).
The two Maxent models produced consistent results and had high predictive power (AUCm1 0.98 and AUCm2 0.99; Figure 2, Supplementary Figure S3). The test omission rate of m1 (subsample approach) is close to the predicted omission, while the omission rate of m2 (bootstrap approach) slightly deviates from the predictive values for cumulative thresholds <0.5 (Figure 2, Supplementary Figure S3).
Hot springs with a temperature between 45 °C and 100 °C and the amount of precipitation in the coldest quarter of the year (i.e., winter snowfall; BIO19) had the highest percentage contribution to the models (58%, 15% in m1, and 56%, 14% in m2, respectively; Figure 2, Supplementary Figure S3, Table S2). Both maps of the two models show very similar patchy distribution patterns of habitat suitability with high values (yellow to red) corresponding to pixels that overlap with several hot-spring sites in the central Hengduan Mountains and in the tributary valleys of the Yarlung River in southern central Tibet (Figure 2, Supplementary Figure S3). These suitable habitats are located in areas with lower winter precipitation than surrounding regions (Figure 3; Supplementary Figure S4). Hot springs outside these ranges are associated with low values of habitat suitability (green), and areas without hot springs are unsuitable (blue).

Habitat Characteristics
We detected T. zhaoermii in 6 out of the 22 surveyed sites, including the type locality in Litang county (Supplementary Figure S1; Table S1). All of these few "positive" sites were located west of the Yangtze River at elevations between 3600 and 4090 m (a.s.l.). Three further sites were considered as "uncertain" with respect to the presence of the species and were excluded from subsequent analysis. At these sites, we observed snakes only from long distances in difficult mountainous terrain and could not verify them morphologically.
Spearman correlation coefficients >0.5 for all pairs of variables measured at the hotspring sites ranged between 0.52 and 1.00. Variables TDS, SAL, and NH3/NH4 were excluded from subsequent analysis due to high pairwise correlation (>0.7) with other parameters (Supplementary Table S3); S2¯ was removed due to invariability. The final dataset contained 14 variables. Forward model selection among these variables did not reveal an important contribution of any hot-spring water characteristics, the distance of a hot spring to the nearest small river or to the superordinate river. The final model comprised only two variables: MRD and TEMP (Supplementary Table S4). The former was nominally insignificant in a single-variable likelihood-ratio test (LRT), while TEMP turned out to be highly significant (p = 8.73 × 10 −05 ; Supplementary Table S4), suggesting an impact of the hot-spring temperature on the presence of T. zhaoermii. Mean hot-spring water temperature at sites where hot-spring snakes occur was about 56 °C (41-70 °C), i.e., ~10° K warmer than hot-springs at which the snakes did not occur; the mean distance from these "presence" sites to a major river was ca. 70 km, in contrast to distances >120 km from sites without snakes (Supplementary Tables S5 and S6).

Discussion
Determining how species are distributed in space is a key issue in ecology and conservation biology and this knowledge is particularly important in endangered species.
Our study provides first information on the potential geographic distribution range of the rare species T. zhaoermii (and T. shangrila) and largely confirms previous data on the occurrence area of T. baileyi. The most important predictors in our niche models were hot springs with temperatures ranging between 45 °C and 100 °C, and winter snowfall. Moreover, results of the fine-scale modelling support the relevance of the hot-spring temperature for the presence or absence of Thermophis. They also indicate an impact of the distance from a hot-spring site to a major river on the distribution pattern of the snakes. Although based on a relatively small data set, our study presents valuable insights into the potential distribution of Thermophis, which can assist in setting conservation priorities for those snakes and should encourage future ecological studies in the genus.

Habitat Characteristics
We detected T. zhaoermii in 6 out of the 22 surveyed sites, including the type locality in Litang county (Supplementary Figure S1; Table S1). All of these few "positive" sites were located west of the Yangtze River at elevations between 3600 and 4090 m (a.s.l.). Three further sites were considered as "uncertain" with respect to the presence of the species and were excluded from subsequent analysis. At these sites, we observed snakes only from long distances in difficult mountainous terrain and could not verify them morphologically.
Spearman correlation coefficients >0.5 for all pairs of variables measured at the hotspring sites ranged between 0.52 and 1.00. Variables TDS, SAL, and NH 3 /NH 4 were excluded from subsequent analysis due to high pairwise correlation (>0.7) with other parameters (Supplementary Table S3); S 2− was removed due to invariability. The final dataset contained 14 variables. Forward model selection among these variables did not reveal an important contribution of any hot-spring water characteristics, the distance of a hot spring to the nearest small river or to the superordinate river. The final model comprised only two variables: MRD and TEMP (Supplementary Table S4). The former was nominally insignificant in a single-variable likelihood-ratio test (LRT), while TEMP turned out to be highly significant (p = 8.73 × 10 −05 ; Supplementary Table S4), suggesting an impact of the hot-spring temperature on the presence of T. zhaoermii. Mean hot-spring water temperature at sites where hot-spring snakes occur was about 56 • C (41-70 • C), i.e., 10 • K warmer than hot-springs at which the snakes did not occur; the mean distance from these "presence" sites to a major river was ca. 70 km, in contrast to distances >120 km from sites without snakes (Supplementary Tables S5 and S6).

Discussion
Determining how species are distributed in space is a key issue in ecology and conservation biology and this knowledge is particularly important in endangered species. Our study provides first information on the potential geographic distribution range of the rare species T. zhaoermii (and T. shangrila) and largely confirms previous data on the occurrence area of T. baileyi. The most important predictors in our niche models were hot springs with temperatures ranging between 45 • C and 100 • C, and winter snowfall. Moreover, results of the fine-scale modelling support the relevance of the hot-spring temperature for the presence or absence of Thermophis. They also indicate an impact of the distance from a hot-spring site to a major river on the distribution pattern of the snakes. Although based on a relatively small data set, our study presents valuable insights into the potential distribution of Thermophis, which can assist in setting conservation priorities for those snakes and should encourage future ecological studies in the genus.
Species of the genus Thermophis are the only snakes with a vertical distribution up to 4900 m a.s.l. [29]. The only snake taxon to come close in the Himalaya is Gloydius himalayanus which has been recorded at 4772 m [45]. Living in such extreme habitats requires specific physical and or behavioral adaptations to cope with the high, cold, and arid environments. Besides molecular adaptive responses to oxidative stress and UV radiation that have been shown in Thermophis [21], the strong association of these snakes to hot-spring habitats has been interpreted as an ecological adaptation that evolved during the geological uplift of the Tibetan Plateau [19,29]. Irrespective of the variable taxonomic placement of Thermophis (Pseudoxenodontinae vs. Dipsadinae; [46] and references therein) and wide range of its estimated divergence time (MRCA Thermophis and xenodontine snakes~10-28 Mya; [47]), ancestral lineages of this genus might have been present in the area of Paleo-Tibet long before the final uplift of the Plateau. With the continuously rising Himalaya-Tibet orogen and the associated cooling of the environment, these ancestors have probably been forced to retreat into geothermal active areas, or, otherwise, would have gone extinct. At least some of the hot-spring sites have served as refuges for the snakes during the Pleistocene [19]. It can be assumed that in the alpine zone hot springs may have become essential for Thermophis to be present to survive during cold periods and to achieve enough activity time that determines the species' fitness [48]. However, not every hot-spring site represents a suitable habitat for Thermophis. The high predictive power of springs with water temperatures between 45 • C and 100 • C in our models suggests a clear preference of the snakes for moderately but not too warm places. We suspect that the geothermal sites are particularly important to the snakes during and soon after hibernation in wintertime, and less relevant throughout the highest activity peak in summer. This may also explain why Thermophis can occasionally be found more distant from a hot-spring site during the middle of the year ( [49]; pers. obs.), especially throughout the mating season [50]. Since hibernation sites at high elevations must protect snakes from colder conditions for longer periods of time, we assume that Thermophis hibernates underground near the geothermal field in rodent burrows, old root systems, or rocky outcrops. Here, they probably prefer thermal sites that are warm enough to protect them from the cold, but not hot enough to achieve their body temperatures necessary for appetite and digestion. Geothermal fields with very hot water temperatures may heat up deep soil layers that then become too warm and unsuitable for hibernating snakes.
The broad-scale influencing climatic factor for the distribution of Thermophis taxa is linked to precipitation during the coldest quarter of the year, with an increase in winter snowfall decreasing the probability of the snake's presence in a hot-spring habitat. Therefore, regions suitable for Thermophis experience significantly lower precipitation during the winter than do areas predicted to be unsuitable for the genus. A similar finding has been reported for Eirenis persicus in western Iran and Turkey [51]. In Thermophis, an intuitive explanation for this pattern is that the advantage of the snowpack by insulating the belowground hibernation site from ambient temperatures [52,53] might be negated in geothermal active locations. Moreover, due to the generally long hibernation season, high-altitude reptiles depend on the relatively short active season to maximize growth, reproduction, and energy storage. A higher winter precipitation, however, can result in the accumulation of a large seasonal snowpack, prolonged cold and, thus, an even shorter activity period for the snakes.
At a fine scale, our results support the hypothesis that the hot-spring temperature, and probably also the access to the hot spring over short distances via a major river, influence the distribution of Thermophis, although these findings have to be treated with caution due to the low number of observations. If true, long riverine pathways may have prevented a suitable hot-spring site being colonized by Thermophis in the past. It is well-known that river-mediated migration can be a key factor for the distribution of species, for example in the European pool frog Rana lessonae [54], and the common wall lizard, Podarcis muralis [55]. In a previous molecular study, we could show that populations of T. baileyi are connected through riverine corridors by the Yarlung River and its tributary system [17]. We suggest that, for hot-spring snakes, the suitability of a geothermal habitat may depend on the trade-off between species dispersal capacity and the prospects for successful settling that could be decreasing with increasing along-valley distance to the hot spring.

Conclusions
Today, the main areas of occurrence for Thermophis cover the southern central region of the Tibetan Plateau along the Yarlung suture zone, and the central range of the Hengduan Mountains (Chola Shan, Shaluli Shan). The fact that only these relatively narrow areas are predicted to be suitable, and host recorded presence localities suggests that the potential distribution of the species might not be much wider than our current records indicate [17]. The distribution map of the genus presented here provides details on hot-spring localities with a high probability of the species' presence, which could guide future surveys and conservation activities. It also indicates that the current distribution of these snakes depends on the amount of winter precipitation, while other climatic predictors seem to be less informative with respect to the distributional patterns of the genus. According to the IUCN Red List T. baileyi is considered as "near threatened", while T. zhaoermii as "endangered"; T. shangrila is not even registered in the IUCN red list. Both T. baileyi and T. zhaoermii are flagged with a decreasing population trend and their presence in any protected areas is unknown [15]. Given the geographic distribution maps in the most recent status assessment for the IUCN Red List, our predicted potential distribution of Thermophis differs substantially from these data. However, despite their exceptional value for biodiversity, no specific conservation action plans exist for the hot-spring snakes. Our results contribute to the knowledge about the distribution of these species and may guide future management of the genus, e.g., by surveying and protecting suitable habitats for the hot-spring snakes.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/d13070325/s1, Text S1: Details on site visits. Table S1: Hot-spring localities surveyed for Thermophis zhaoermii. Water parameters were measured at hot-spring sites 1-16 (see Figure S1 for geographic reference, below); localities 17-22 were surveyed for snakes only. Table S2: Estimates of relative contributions of the environmental variables to the Maxent model based on subsample (left) and bootstrap algorithm (right). Table S3: Spearman correlation coefficients (SCE) of microhabitat characteristics. The following variables showed pairwise correlation > |0.5|. Table S4: Analysis of deviance table (generalized linear model). Table S5: Summary statistics of variables in the model for "absence" localities of Thermophis zhaoermii. Table S6: Summary statistics of variables for "presence" localities of Thermophis zhaoermii. Figure S1: Known records of Thermophis zhaoermii and T. shangrila and sites surveyed for Thermophis. Details to the locality numbers are listed in Table S1. Figure S2: Climate heterogeneity in the distribution area of Thermophis species. Figure S3: Diagnostic plots for the Maxent model of Thermophis based on bootstrap algorithm and 100 replicates. Figure S4