Historical and Ecological Factors Affecting Regional Patterns of Endemism and Species Richness: The Case of Squamates in China

1. Introduction

Biogeography is closely tied to both ecology and phylogenetic biology and its main areas of interest are ecological biogeography, i.e. the study of factors influencing the present distribution, and historical biogeography, i.e. the study of causes that have operated in the past [1]. Ecological and historical biogeography therefore applies different concepts in order to explain the distribution of organisms. The former deals with functional groups of species and environmental constraints, whereas the latter focuses on taxonomic groups and historical biogeographical events [2]. Recently, the division between historical and ecological biogeography has been considered as an obstacle to the progress of biogeography and some authors have stressed the benefits of integrating these two points of view [1,3-4]. In this context, the present work attempts at integrating these two approaches to explore the regional patterns of endemism and species richness of the squamates(lizards and snakes) in China.

As summarized by Meng et al. [5], China has a relatively independent geological history. Six primary plates are involved in its tectonic history, namely the North China, Tarim, Yangtze, Cathaysian plates, and parts of the Siberia and Gondwana plates [6-7]. In the Pleistocene, seven collisions and integrations of these plates had united the ancient Siberian and European plates. With approximately 200 tectonic-facies, China has a complex topography, including towering mountains, basins of various sizes, undulating plateaus and hills, and fertile plains. As noted by Meng et al. [5], they are assigned to four different terraces in general. The highest terrace, the Qinghai-Tibetan Plateau – the ‘roof of the world’ – has an average elevation of over 4000 m. Tremendous differences in latitude (a span of more than 50°latitude), longitude (a span of more than 60°longitude), and altitude (a span of more than 88km) create the conditions for extremely diverse climate and a highly astounding heterogeneous landscape. Thus, a high degree of species richness and endemism in China as might have been expected in terms of squamates diversity. It is conservatively estimated that today there are about 422 squamates species in this country, belonging to 106 genera in 17 families and two suborders (see Appendix S1 in Supplementary Material).

It has long been recognized that geological complexity and history usually have a profound influence on the distributions of living organisms. China animal geography division was firstly put forward in 1959 [8]. Now based on the distributions of vertebrates, mainly mammals and birds, China was divided into seven major biogeographical regions and fifty-four provinces [9]. The Palaearctic realm includes North-eastern, Northern, Inner Mongolia-Xinjiang, and Qinghai-Tibetan China biogeographical regions, while the Oriental realm was divided into South-western, Central, and Southern biogeographical regions. The study on biogeographical divisions of China was both intriguing and challenging due to the complex topography and historical processes [9]. The zoogeographic division of China was explored into 30 units [10]. However, China’s territory was divided into 124 basic units on the basis of comprehensive natural factors (including altitude, landform, climate, vegetation, water system, farm belt, and so on) by cluster analysis [11]. Recently, based on the distributional patterns of spiders, seven major biogeographical regions in China were investigated by PAE, not corresponding to any previous studies [5]. These studies played an important role in biodiversity conservation, planning and management in China. However, these divisions need further investigation due to primarily on basis of mammals and birds which are highly adapted to environment diversity and higher locomotion [12]. Thus, the ectothermic animals may act as a better indicator to determinate zoogeographical division than mammals and birds [12]. Although only a few studies on biogeographical patterns have been identified and emphasized their biogeographical complexity, data are still not enough to evaluate and compare directly for other taxa, such as squamates with limited dispersal abilities. It is necessary to compare patterns based on different groups of organisms to better understand their biogeography and infer a general pattern [13].

2. Compared biogeographical patterns of Chinese squamates based on parsimony analysis of endemicity (PAE) at different natural area units

2.1. Introduction

One of important goals of historical biogeography is to investigate convergent biogeographical patterns relying on different taxa [14]. These may assist in identifying priority areas or hotspots for biodiversity conservation, particularly today the issue of global concern biodiversity loss. A historical biogeographical method, parsimony analysis of endemicity(PAE), firstly proposed by Rosen [15] and further elaborated by Morrone [16], provided an insight to generate area cladograms to make inferences on historical patterns. PAE was originally used the most parsimony algorithm to reconstruct relations among sampling localities [15], then previously delimited areas [17] and quadrats [16]. Analogous to cladistic methods in phylogenetic analyses, PAE classifies areas (cf. taxa in cladistics) on the basis of the shared presence of taxa (cf. characters in cladistics) [15]. Using PAE, biogeographical studies can investigate biotic similarities between different geographical regions, give static or ecological interpretations [18], and estimate historical hierarchical congruence in target localities or geographical regions [19]. Although there is ongoing debate about the value of PAE [20-21], the PAE has been widely used in many biogeographical studies in recent years [5,22-23].

PAE has been applied to establish relationships among different areas units, for example localities, quadrats, areas of endemism, continents, islands, and so on [13]. The ideal organisms for using PAE are those with limited dispersal abilities and speciation in vicariant events [13,19]. Although the same size and shape quadrats are not required in PAE and do not affect the analyses [24], the best PAE results were obtained with natural areas instead of quadrats [13]. Here, we used PAE to compare squamates biogeographical patterns at different natural area units with previously delimited biogeographic patterns obtained by Meng and Murphy [5] (Figure1), Zhang [9](Figure 2) and Xie et al. [11](Figure 3), in order to find areas of congruent distributional patterns in China.

2.2. Materials and methods

2.2.1. Study area and operational geographical units

For a comparable study area with previously provided optimal results in detecting the biogeographical patterns, we assumed the same operational geographic units (OGUs) as suggested by Meng and Murphy [5] (Figure 1), Zhang [9] (Figure 2) and Xie et al. [11] (Figure 3) defined, respectively. According to very similar climatic, geological characteristics, topographical characteristics and natural barriers to dispersal for neighbouring quadrats, 28 biological and physiographical similar areas were combined and divided as OGUs in PAE [5]. Zoogeographical divisions of China were classified into 54 biogeographical provinces based on the distribution of mainly mammals and birds [9]. The whole China was divided into 124 basic units using comprehensive natural factors, such as altitude, landform, climate, vegetation, water system, farm belt, and so on [11]. More details about the OGUS see these studies.

2.2.2. Data sources and dataset

The Squamates species catalogue and distributions were compiled from the most recent and comprehensive references, including specimens, exhaustive field surveys, monograph, published literature and expert reviews. All date Information from herbarium specimens was mainly obtained through Chengdu Institute of Biology, Chinese Academy of Sciences (CIB/CAS). Additional information was also obtained from HerpNET and a variety of published sources, for example Fauna Sinica Reptilia, Vol.2 Squamates Lacertilia [25], Herpetology of China [26], Fauna Sinica Reptilia, Vol.3 Squamates Serpentes [27], Snakes of China (2 Volume set) [28], Zootaxa, Zoological Journal of the Linnean Society,Bonn Zoological Bulletin, Asian Herpetological Research, Zoological Research, Acta Herpetologica Information, Acta Herpetologica Sinica, Acta Zootaxonomica Sinica, Sichuan Journal of Zoology, Chinese Journal of Zoology, and so on. For the taxonomic revisions and updated information recently, we follow the taxonomy of Fauna Sinica Reptilia [26-27], Snakes of China [28] and the reptile database (http://www.reptile-database.org/). To reduce the sources of potential errors records, uncertain distributions and identifications, such as isolated distribution records, we checked again and removed erroneous records. If possible, we tried to overcome these difficulties by adding as many taxa as possible to the analyses. The distributional data with longitudes and latitudes were reviewed from museum databases or literature, and otherwise, we compiled coordinates from Google Earth. Finally, an initial data set of 422 Squamates with species names (see Appendix S1 in Supplementary Material) represented by a total of 50,749 records, of these, 61.0% records were lizards, while the remaining records were snakes.

2.2.3. Parsimony Analysis of Endemicity (PAE)

We used PAE to identify biogeographical patterns and followed the procedure modified by Morrone [16]. A taxon/area matrix for the basic PAE data set was built in which the absence of a species in an area was coded as ‘0’ and presence as ‘1’. A hypothetical area coded ‘0’ as an outgroup to provide a root for the final cladogram [13, 16]. We removed species that occurred in all areas, as well as those in only one area due to phylogenetically uninformative autapomorphies [15]. The taxon/area matrix were imported into PAUP * [29] to find the most parsimonious cladograms with a heuristic search of 1000 replicates and random sequence additions. We estimated relative support for each branch using bootstrapping, with 100 replicates and tree bisection–reconnection (TBR) swapping. All characters were weighted equally and a 50% majority consensus tree of the equally parsimonious trees was generated.

2.3. Results

2.3.1. Biological areas 28

The bootstrap 50% majority-rule consensus tree using 244”characters” in the analysis was shown in Figure 1(b). Our results were substantially different from the results of Meng et al [5]. The parsimony analysis obtained twelve most-parsimonious trees [Tree length = 676; Consistency index (CI) = 0.3609; Homoplasy index (HI) = 0.6391; Retention index (RI) = 0.6415]. Although two clear groups, D and E, were identified by the analysis, they did not support either the delimitation between palaearctic and oriental realms in China or the geographical barrier within the confines of the Qinling Mountains and Huai River, and east of the Hengduan Mountains [5]. The bootstrap value of group D (64%) was almost the same as that of group E (61%). The group D includes area 21 and area 23, which corresponded to the Tibetan Plateau and Taiwan Island in Map (Figure 1).

For Group E, it was very complex and ten clades were discovered. Subclade E1 as Sothern region (bootstrap 88%) corresponded to the clade C in [11], but major regions such as Central region (C2), Eastern Southern region (C3), Western Southern region (C4), and Central Southern region (C5) in [5] cannot be recognized. In subclade E1, A8 mainly representing Jiangsu was basal, followed by separation of A1 and the remaining areas in subclade E1. Areas A22, A24, A25, A26 and A27 formed a clade with a moderate support of 78%, corresponding to Southern China of zoogeographical division of [9] excluding Taiwan Island. Subclade E2 had a bootstrap support of 98%, including two sub-clades, A12 and A5 + A28. It corresponded to northeastern Tsaidam basin, Loess plateau and Alashan Plateau. Subclade E3 (A10 + A19) corresponded to the steppe and desert of north-western China, including Xinjiang and Inner Mongolia. Its bootstrap support value was 72%. Subclade E4 (A14 + A18) consists of the Xiao Xingan and Changbai mountains, with a bootstrap value of 86%. The Qinghai-Tibetan Plateau included A11+A21, however, subclade E6 (A11) was separated from A21 in the results with a bootstrap support value of 64%. The remaining subclades showed ambiguous.

Figure 1.

(a) The division of 28 areas (A01–A28) in China used in the PAE analysis. (b) 50 % majority consensus tree generated by PAE analysis. The terminals correspond to the OGUs shown in (a). Bootstrap values were shown above the branches

2.3.2. Biological areas 54

The 1507 most parsimonious cladograms of 930 steps with CI of 0.3172, HI of 0.6828, and RI of 0.6340 were found. The 50% majority-rule consensus tree using 295”characters” (Figure 2, b) showed a basal polytomy, but several clades emerged. Clade F consisted of areas A12, A13 and A14 with a bootstrap value of 57%, which corresponded to the Loess plateau subregion. Clade G had a weakly supported group containing areas A18, A19, A20, A21, A23 and A24 with a bootstrap proportion of 57%, which corresponded to Western desert subregion and Tian Shan mountains subregion of Mognolia-Xijiang China. Clade H had a well-support of 95%, including two sub-clades, A30 + A36 and A37, in response with the Himalaya mountains subregion of Southwestern China. Areas A40, A44, A45, A46, and A47 formed a clade I, including Eastern hills and plains subregion and Western mountains and plateau subregion of Central China, and coastal subregions of Guangdong and Fujian provinces of Southern China; its bootstrap value was 57%. Clade J consisted of areas A48 and A49 with a strongly support value of 100%, corresponding to Southern Yunnan mountains subregion of Southern China. Clade K corresponded to Hainan Island subregion of Southern China, consisting of areas A50 and A51 with 75% bootstrap value. Clade L contained Taiwan Island areas with similar Squamates faunas, areas A52 and A53; its well-supported value was 100%.

Figure 2.

(a) The division of 54 areas (A01–A54) in China used in the PAE analysis. (b) 50 % majority consensus tree generated by PAE analysis. The terminals correspond to the OGUs shown in (a). Bootstrap values were shown above the branches

2.3.3. Biological areas 124

For OGUs of biological areas 124, the 1000 most parsimonious trees were obtained. The 50%-majority consensus tree using 294 characters (tree length = 1407, CI = 0.2090, RI = 0.7910) of the 1000 trees was shown in Figure 2 (b). Like the Figure 3, a basal polytomy were identified where only one of the branches contained a dicotomy in Figure 3. Clade M is moderately supported by Areas A10 and A14 with 88% bootstrap value. It corresponded to Yinshan mountains hill and Ordos Plateau. Areas (A21, ((A22, A23), A17)) formed the clade N with a week support value of 62%, corresponding to (Mountain front mesa, ((Liaodong Peninsula, Lower Liaohe Plain), Changbai mountains)) in Eastern Northeastern China. Clade O consisted of areas A29 and A32, including Southeast Shanxi plain and Western Henan mountains. Its bootstrap support proportion was 60%. Clade P contained areas (A101, (A37, A99)) with 68% bootstrap value, corresponding to the range of Qilian Mountain, Central Gansu incisive hill and Upper Yellow River incisive mountains. Areas A46 (Qinling mountains) and A47 (Daba–Micang mountains) in Central China formed clade Q with a moderate support value of 81%. Clade R was compose of Areas A58 and A63, corresponding to Xiangjiang valley hill and Honghe catchment montane basin; its support proportion was 80%. Clade S included Areas A61, A62 and A63 with a weekly support value of 51%, corresponding to Dalou Mountain mid-land valley, Miaoling hilly plain and Wujiang and Nanpanjiang catchments mid-land valley. Clade T was the largest clade including Areas A67, A68, A73, A74, A75 and A87 with 60% support value. It corresponded to Southeast Yunnan low-heat plateau, South central Yunnan low-heat valley, West Yunnan montane plain, Southwest Yunnan Plateau wide valley, South Yunnan wide valley and Salween and Lancang Rivers parallel valley. Areas A76 (West Guangdong and south Guangxi coastal mesa plain) and A79 (South Hainan montane hill) formed Clade U with 64% support value. Clade V had the greatest bootstrap value (98%), including areas A80, A81 and A83. It corresponded to Northwest subtropical hilly plain, Central subtropical mountain and East tropical coast in Taiwan Island. Clade W consisted of areas A92, A93, A94 and A95 with a bootstrap support of 73%, corresponding to areas in Himalayas mountains including Kangrigebu south wing mountains, Himalayas south wing mountains, Salween and Lancang Rivers incisive mountains and Brahmaputra Great Turn and upper Salween incisive mountains. Similarly, clade X consisted of areas A96, and A97, corresponding to areas in Himalayas mountains including Brahmaputra valley mountains and Himalayas central mountains; its bootstrap support value was 61%. Clade Y was compose of areas A104 and A105 with a week bootstrap support proportion of 59%. It corresponded to areas in Qinghai–Tibetan Plateau including South Qiangtang Plateau mountains and North Tibet plateau northwestern lake basin mountains. Clade Z included areas (A111, (A115, (A120, A122))), which corresponded to West Hexi Corridor, East Tianshan mountains, Central Tianshan mountains, Tarim Basin. Its bootstrap value was 57%. Areas A118 (Junggar Basin), A121 (Ili Valley), and A119 (Emin Valley) formed Clade AZ with a well support proportion of 97%.

2.4. Conclusion

Past major geological events have played important roles in shaping the biogeographic distribution of extant organisms. PAE originally aimed to find areas of congruent distributional patterns, and the best PAE results were obtained with natural areas (e.g. biogeographical provinces, ecoregions) instead of quadrats by increasing the absolute and relative numbers of synapomorphies [13]. If we compare the results of areas 28, 54 and 124, there seems to exist a trend to result in poor resolution of the resultant area cladograms as the size of OGU decreases. Although the samples within the regions are not uniform at different areas unit, in our area cladograms, biogeographic patterns of squamates distributions appear to have a hierarchical structure and general patterns for areas 28. Based on the above comprehensive squamates distributions patterns at different natural area units, seven major congruent biogeographic regions can be identified in China: Eastern Northern region, Tibetan Plateau region, Xinjiang and Inner Mongolia region, Loess plateau and Alashan Plateau region, Taiwan Island region and Southern region. However, there existed several unresolved areas relationships to one another, the uncertain position of certain areas or incongruence in the ‘character’ distributions. There may be basal problems with data themselves such as squamates distribution incomplete known. Sure, extinction, long-distance dispersal, isolation, or other undetected historical patterns may lead to incongruence in the distributions patterns [30].

Figure 3.

(a) The division of 124 areas (A01–A124) in China used in the PAE analysis. (b) 50 % majority consensus tree generated by PAE analysis. The terminals correspond to the OGUs shown in (a). Bootstrap values were shown above the branches

It is surely not coincidental that different organisms may share a general distribution pattern. A PAE area cladogram might contain areas related by shared ecologies or similar historical events (biotic divergence and isolation) [15]. Historical hypothesis, in other words, evolutionary history has recently been considered to be a driving force determining squamates regional species pools’ differences in China. Geological complexity and history usually have a profound influence on the distributions of living organisms in China [5], though most of China has never been covered by ice sheets. Squamates species are followed this rule. For example, squamates were similar to spiders distributional patterns broadly corresponded to geological provinces in China, such as Southern China geological province versus Southern regions, The Laurentian/Cathaysian Southern, South-western Margin geological province and Tibetan geological provinc vs Tibetan Plateau biogeographical region [9]. Furthermore, biotic and abiotic conditions are also important factors in determinative of the distribution of squamates species. Some correlations between species richness and reproductive modes with geography and ecological conditions have been reported [31-32]. The interpretations of those relationships have postulated that contemporary factors are the main regulatory force of the distribution of squamates taxa in China [33-35]. Thus, both history and ecology may well be inseparable and have a profound impact on not only the diversity of Squamates taxa but also their biogeographic patterns.

3. Distribution patterns in species diversity of lizards in China and their relationships to ecological factors

3.1. Introduction

Determining the causes of the great biodiversity variation across Earth has long been a major challenge for ecologists and biogeographers [36], ever since biotic diversity contrast between equatorial and polar latitudes was discovered two centuries ago [37]. Among the considerable number of hypotheses that aim to explain species richness patterns [36, 38], many ecological (environmental) hypotheses have been widely discussed and accepted [39].

Three alternative variants of ecological hypothesis, the species-energy, contemporary climate and habitat heterogeneity hypotheses, have received a great deal of attention as the primary determinants of species richness [39-43]. The species-energy hypothesis includes at least two versions, the ambient energy and productive energy hypotheses [38]. The ambient energy hypothesis, widely indicated by temperature or allied measures, argues that species richness was influenced by energy inputs into an area that affects the physiological tolerance of organisms [40, 44]. The productive energy hypothesis claims that animal species richness is limited by energy via food webs rather than by physiological requirements. The energy and water availability (i.e., energy–water dynamics) limits the total available plant productivity, which ultimately moves up the food chains [40, 45-47]. The contemporary climate hypothesis states that species richness correlates with contemporary climate conditions, and putative causal mechanisms are in terms of environmental stability, variability, favorability and harshness [40, 47-48]. The habitat heterogeneity hypothesis is measured either as the number of habitat types or the topographic relief (range in elevation) presented within an area [42, 50]. It assumes that high species richness is found in physically or biologically complex habitats, through higher speciation rates and providing more ecological niches [40-42].

However, the knowledge of the determinants of reptile richness remains insufficiently documented among terrestrial vertebrates [39, 51-82]. It is urgent to understand the drivers of reptile richness patterns due to global warming impact on species distribution and abundance [36, 52-53]. Lizards belong to Reptila and are good model systems to test these alternative hypotheses. Because their taxonomy is well resolved and distributional data are quite thorough. They are ectothermic and sensitive to environmental variables. In this study, we examine the correlation between lizard species richness and various environmental factors across China. Our objectives include (1) mapping distributions of Chinese lizards and describing any patterns, and (2) testing various ecological factors in determining species richness patterns.

3.2. Materials and methods

3.2.1. Data collection

We collected locality data for lizard species which occur in China from a variety of sources as above-mentioned. We excluded coastal grid cells with less than 96% land cover and all islands from the analysis in order to remove the effects of insularity. Finally, we built a database of 151 lizard with species names (see Appendix S1 in Supplementary Material) represented by a total of 3,391 records for unique point localities, with a range of 2–288 (mean = 22.5, standard deviation = 38.3).

3.2.2. Ecological niche modeling and species richness

For each of the 151 species, we used the Genetic Algorithm for Rule-set Prediction (GARP) (for free download see: http://www.nhm.ku.edu/desktopgarp/) [54] for reconstructing species distribution maps. GARP uses an evolutionary computing genetic algorithm to search iteratively for non-random correlations between species presence and environmental variables for localities using several different types of rules (i.e., atomic rules, range rules, negated range rules, logistic regression rules), and then creates ecological niche models for each species’ predicted distribution, as contrasted with environmental characteristics across the overall study area [54]. GARP was found that it did not tend to be more sensitive to sampling bias than Maxent, and GARP is a very useful technique to estimate richness and composition of unsampled areas and have been tested to correctly predict the most of the species’ distributional potential [48, 55-58], for example in applications to invasive species [57-58], tree species [59-60], squamate species [48, 61], and so on.

We included a total of eighteen environmental variables in the model. Variables for details, descriptions, and files for download are described in the following text. We set several optimization parameters while running the software following [48]. The parameters included: 20 runs, 0.001 convergence limit, and 1,000 maximum interactions; rule types: atomic, range, negated range, and logistic regression; best subset active, 5% omission error, 40% commission error, and 67% of points for training; omission measure = extrinsic, and omission threshold = hard; 10 models under hard omission threshold.

The estimation output of DesktopGarp produced in Arc/Info grid maps with ‘zeros’, where the species were not predicted to occur, and ‘ones’, where the species were predicted to occur. The area covered by the coincidence of at least seven out of the 10 models in the best subset selection (optimum models considering omission/commission relationships [62]) were used as the predicted distribution of each species. By doing so and by setting the commission error to 40%, this approach added a component of conservatism in predicting distribution by GARP, which might otherwise extrapolate too much and predict areas that are too far from where the species have previously been collected [48]. After generating such maps using the same criteria for all 151 species, we used ARCGIS software to overlay all species prediction maps into a composite map. This final map was used to create a girded of species richness map at a resolutions of 100 km (approximately equivalent to 1°at the equator) on an Albers Equal-Area Conic projection. Consequently, we used the occurrences of 151 lizard species within 827 grid cells to calculate species richness, summing the value of overlaid corresponding grid cells.

3.2.3. Environmental data

We used eighteen environmental variables. We selected these variables based on previous studies and the four associated hypotheses [40-43, 49]. All environmental variables for assessing hypothesized explanations of species richness were re-projected and re-sampled to the same equal-area cell as the species richness data in ARCGIS. The hypotheses and their related variables are:

1. Ambient energy—five variables are associated this hypothesis within each cell, including: mean annual potential evapotranspiration (PET) ([63], 30'resolution, available at http://www.grid.unep.ch/data/grid/gnv183.html); mean annual highest temperature (HT), and mean annual lowest temperature (LT) (data from 1961 to 1990 with 1 km² resolution, available at http://www.data.ac.cn/ index.asp); mean annual sum of effective temperature (≥0℃) (SET0) and mean annual sum of effective temperature (≥10℃) (SET10) (data from 1981 to 1996 with 500 m² resolution, available at http://www.geodata.cn/Portal).

2. Productive energy—three variables are used to account for productive energy hypothesis, including: mean annual remotely sensed Normalized Difference Vegetation Index (NDVI), obtained from Advanced Very High-Resolution Radiometer (AVHRR) record of monthly changes in the photosynthetic activity of terrestrial vegetation (data from 1998 to 2008 with 1 km² resolution, Data source: Environment and Ecology Scientific Data Center of western China, National Natural Science Foundation of China, available at http://westdc.westgis.ac.cn), mean annual actual evapotranspiration (AET)( [60], 30' resolution, available at http://www.grid.unep.ch/data/grid/gnv183.html), and mean annual solar radiation (RAD) (data from 1950 to 1980 with 1 km² resolution, available at http://www.geodata.cn/Portal).

3. Contemporary climate hypothesis—eight variables are associated with this hypothesis within each cell, including: mean annual temperature (AT) (data from 1961 to 1990 with 1 km² resolution, available at http://www.data.ac.cn/index.asp); mean annual sunshine (SUN) (percent of daylength), mean annual diurnal temperature range (DTR) and mean annual frost-day frequency (FF) (data from 1961 to 1990 with 10' resolution [64]); and mean annual wind speed (WIND) (data from 1981 to 1996 with 500 m2 resolution, available at http://www.geodata.cn/Portal); mean annual precipitation (PRE) (data from 1961 to 1990 with 1 km² resolution, available at http://www.data.ac.cn/index.asp), mean annual wet-day frequency (WET) (number days with >0.1 mm precipitation per month) and mean annual relative humidity (REH) (data from 1961 to 1990 with 10' resolution [64]).

4. Habitat heterogeneity—the count of 300 m elevation range within each quadrat (ELE) (HYDRO1 k data set for Asia, 1 km² resolution, available at http://eros.usgs.gov/) and the number of vegetation classes (VEG) (1 km² resolution, Data source: Environment and Ecology Scientific Data Center of western China, National Natural Science Foundation of China, available at http://westdc.westgis.ac.cn) as indicators of habitat heterogeneity.

3.2.4. Statistical analyses

In order to examine the potential predictors of lizard richness patterns in China, we first tested the relationship between lizard richness and environmental variables using a multiple regression analysis. We did not use all environmental variables employed to run GARP, because including many highly correlated variables in a multiple regression creates several theoretical and statistical problems, especially estimating partial regression coefficients [65]. We selected variables previously identified as affecting species richness and were not highly correlated (r<0.80) and there were one variable represented each hypothesis at least.

We used the eigenvector-based filtering, or spatial eigenvector mapping (SEVM) obtained by Principal Coordinates Neighbour Matrices (PCNM) to account for spatial autocorrelation [66]. Spatial autocorrelation is a potential problem when work with large-scale ecological data and explanatory variables [67-68]. Failure to account for spatial autocorrelation could result in inflating Type I error because model fitting may generate artificially narrow standard errors due to the lack of independence among residuals [67-68]. A truncation distance of 102.33 km, calculated in SAM—Spatial Analysis in Macroecology [69], was used to create the spatial filters. Eigenvector filters were chosen when their influence on species richness was both statistically significant (P<0.05) and had sufficient explanatory power (r²>0.02). We selected eigenvector filters in an iterative process, by minimizing both the spatial autocorrelation among residuals and the number of filters used in regression. Moran’s I coefficient were used to examine the model residuals of the spatial autocorrelation in reducing spatial autocorrelation [70]. These filters were then used as candidate predictor variables, together with other environmental predictors formed in the full model. In this way, the effects of environmental predictors are evaluated as partial effects, taking spatial factors into account explicitly [69]. The total explanatory power, r² values, was divided into three parts: a part explained by space, a part explained by environmental variables, and a part of shared explained variance.

To test which hypothesis best explains variation in lizard richness in China, we conducted separate regressions to fit each of the hypothesis, with an addition of mixed models using all variables associated with each hypothesis. The sample-size-corrected Akaike information criterion (AICc) was used to evaluate the goodness of model fit. The model with the lowest AICc score was considered the most parsimonious, therefore optimizing the tradeoff between bias and precision in model construction [71]. The difference between any candidate models and the best model (∆AICc) was used to evaluate the relative model fit when their AICc scores were close. The larger the ∆AICc, the less possible is the fitted model as being the best approximating model in the given models set. In general, Models having ∆AICc ≤ 2 have substantial support (evidence), those in which 4 ≤ ∆AICc ≤7 have considerably less support, and models having ∆AICc≥10 have essentially no support [72]. Model-averaging of estimates using Akaike weights (wi) was used to confront model selection uncertainty [73].

Finally, to make a comparison between the actual available data and the lizard distributions predicted by ecological niche modeling, we mapped species locality points’ data and calculated species richness at 100 km resolution (Figure 4a). This method allowed us to check whether a spatial sampling bias was shown in the final modeling map (i.e., areas that have more species collected coincide with the areas the model indicated as higher species richness) [48].

Figure 4.

Number of species per grid cell based on (a) the raw data of museum collections, and (b) the ecological niche modeling of 151 species. The grid corresponds to the approximate area of China and the area of each cell is 100 km². Blank cells have no specimen based on the major collections [34].

3.3. Results

3.3.1. Patterns of species richness

Figure 4(b) shows the distribution map which summed all 151 lizard species richness in China. Lizard species richness varied between 1 and 38 species per cell (mean: 13 ± 8 SD) and displayed a consistent pattern that species number increased from higher latitudes to lower ones, and from west to east. The Highest species richness occurs in the Oriental Realm tropics, around the border between southern China and southwestern China, and around the Nanling Mountains, at the border between southern China and central China. Other areas with relative high richness included southwestern China, southern China, northwestern and eastern central China of the Oriental realm (Figure 4b). The raw data map shows a slight sampling bias to the northwest central China and northeast southwestern China (Figure 4a), where the largest herpetological museum CIB/CAS in China are located. However, the niche modeling results are not highly influenced by this bias, since areas corresponding to the highest species richness in China do not overlay completely with the pattern. Furthermore, high richness areas were found by the niche modeling in the relatively poor sampling regions, such as the northwest Mongolia-Xinjiang China, west-south Qinghai-Tibet China areas (Figure 4b).

3.3.2. Species richness and environmental variables

The environmental models for the multiple regression analysis using SEVM with adding eigenvector spatial filters (PCNM), were sufficient to reduce autocorrelation in the residuals (filters data not shown). A spatial correlogram based on Moran’s I index was used to evaluate the pattern of spatial autocorrelation in the residuals of the regression. The multiple regression model explained a total of 80.1% variance of lizard richness in China (r² = 0.801; F = 203.47; P<0.001). The total explanatory power explained by predictors alone was 17.0%, and explained by space alone was 5.4%, and the shared explained variance was 57.7%.

Based on the multiple regression analysis taking PCNM spatial filters into account, annual frost-day frequency (FF), elevation range (ELE), the number of vegetation classes (VEG), and wet-day frequency (WET) were the best predictors of species richness. FF was negatively correlated with lizard richness, while ELE, VEG, and WET, respectively, were positively correlated with lizard species richness. Based on model selection approach, the model with the lowest AICc value was the mixed model, which contained all variables related to the different hypotheses. It had an Akaike weight of 1.00. Other models had high ∆AICc values (>10, [69]) and low values of Akaike weights.

3.4. Conclusion

Our results indicate that mechanisms related to different ecological hypotheses might work together to account for lizard richness in China. It is important to consider the influence that environmental factors may have on shaping richness patterns [73]. The frost-day frequency, elevation range, vegetation and wet-day frequency were the most important environmental variable predicting lizard species richness in China. The current alternative hypotheses are not mutually exclusive and may work together and best explain patterns of lizard species richness in China. Based on results of the model selection, our conclusion is in concordant with several previous studies that multiple hypotheses may best account for species richness patterns [48,66,74-75]. Clearly, a variety of factors works synergistically to determine species richness patterns. Our results indicate significant conservation implications, and habitat heterogeneity would be taken into account as an assessment of the threat to endemism from habitat loss in the future investigation. Lizards in China might have experienced large radiations and adapted to dramatic climatic fluctuations after the uplifting of the Tibetan Plateau in Pleistocene. For future studies, it is important to test species richness distribution in Asia at different spatial extent and sample resolution [76], as well as to explore other factors known to affect species richness, such as historical factors and biotic interactions (e.g., competition, predation, and parasitism).

4. Climate and history explain the origin and diversification for toad-headed agamas (genus Phrynocephalus) and racerunner lizards (genus Eremias)

Dramatic geologic events and climatic shifts are often considered to have significantly influenced the diversification and distributions of organisms. Central Asia has a long history of aridity, with the onset of desertification starting at least 22 million years ago (Ma) [77]. It is believed that the Miocene retreat of the Paratethys Sea, an epicontinental sea stretching over Eurasia 30 Ma, and the uplift of the Tibetan Plateau during the Oligocene–Miocene, played a major role in the shift of the Central Asian climate from oceanic to continental [78], leading to increasing levels of aridity. These climatic changes intensified in the Late Miocene and into the Plio-Pleistocene, as part of the global deterioration of the Cenozoic climates [77]. The Quaternary glaciations in the Qinghai-Tibetan Plateau (QTP) and the bordering mountains were the consequence of a combination between climate and local tectonic uplift. In particular, the Kunlun-Huanghe and Gonghe tectonic uplifts have played very important roles in triggering glaciations in high Asia [ [79-81]. The dramatic geological and climatic histories on the Plateau during the Quaternary have a remarkable influence on regional and adjacent biogeographic patterns [82].

4.1. Toad-headed agamas (Phrynocephalus)

The genus Phrynocephalus, i.e., toad-headed agamas, are a modern group of lizards including about ~37 distinct species [83, but see 84]. They are typically grouped together with Asian Rock Agamas (Laudakia), Bufoniceps, and Trapelus in the subfamily Agaminae. Toad-headed agamas inhabit arid regions from northwestern China to the western side of the Caspian Sea, across the Tibetan Plateau, and southwest Asia to the Arabian Peninsula. They constitute one of the major components of the central Asian desert fauna and are highly adapted to sand dunes and desert environments. The reproductive biology of Phrynocephalus is notable in that there exist two reproductive modes: viviparity and oviparity. All six viviparous species are endemic to China and mainly restricted to high elevations in the Tibetan Plateau. Despite the progress that has been made in recent years for the systematics of certain Phrynocephalus groups, the large-scale pattern of their evolution in time and space remains open. Here, we review the origin and diversification of toad-headed agamas with their relationships to geological and climatic changes.

Until recently, our knowledge of the phylogenetic relationships and historical biogeography of Phrynocephalus was mainly based on anatomical studies [85-86]. Researchers usually agree that the broad distribution of Phrynocephalus may be explained by different hypotheses. Ananjeva and Tuniyev [87] speculated that there were two original centers for the Phrynocephalus species: Central Asia (a northern Tethys origin) and Middle Asia (a southern Tethys origin). They also inferred that diversification was related to geological events such as alpine orogenesis causing the isolation of valleys and basins and changes in the direction of river courses. Through a phylogenetic analysis of allozymes, Macey et al. [88] suggested that Phrynocephalus represents an old radiation that has been evolving in response to the Indian collision with Eurasia 35 Ma, and a clade in the former Soviet Union may have diverged as a result of internal Eurasian block movements in Afghanistan caused by the indenting Indian continent. Later, Arnold [86] proposed that Phrynocephalus originated in the Arabia-NW India area rather than in Central Asia, and achieved its current distribution by dispersal. However, Wang and Macey [89] initially considered the origin of the viviparous species group as a result of a vicariance event associated with the uplifting of the Tibetan Plateau, which was subsequently supported by Zeng et al. [90] and Pang et al. [91]. Meanwhile, Pang et al. [91] using mitochondrial DNA data, partial sequences of the mitochondrial genes (12S rRNA, 16S rRNA, cyt b, ND4-tRNA^Leu), corroborated the monophyly of the viviparous group for the first time, and hypothesized that the oviparous group achieved their current distribution by dispersal. Pang et al. [91]. also provided the first time estimate for the origin of Chinese Phrynocephalus at 6.2–3.1 Ma, using average evolutionary rates of rate-constant genes to apply a clock.

Most recent analyses use relaxed-clock methods, which allow evolutionary rates to vary among genes and lineages. Guo and Wang [92] re-analyzed the data of Pang et al. [91] by employing partition-specific modeling in a combined DNA analysis to clarify existing gaps in the phylogeny of Chinese Phrynocephalus. Using this phylogenetic framework, they inferred the genus’ historical biogeography by using weighted ancestral-area analysis and dispersal-vicariance analysis in combination with a Bayesian relaxed molecular-clock approach. They drew three major conclusions: (i) the uplift of the QTP played a fundamental role in the diversification of viviparous Phrynocepalus; (ii) an evolutionary scenario combining aspects of vicariance and dispersal is necessary to explain the distribution of Phrynocephalus; (iii) Chinese Phrynocephalus originated at the Middle-Late Miocene boundary (12.78 Ma, 95% CI 17.61–8.25 Ma), and diversified from Late Miocene to Pleistocene from a center of origin in Central Asia, Tarim Basin, and Juggar Basin temperate desert.

A more recent study [93] using a ~1200 bp region of mitochondrial DNA (ND2-tRNA^Tyr) and a ~1200 bp nuclear gene (RAG-1) and samples across Central Asian agamids results in two well supported conclusions: (i) the onset of aridification in Central Asia during the Late Oligocene, resulting from the retreat of the Paratethys Sea and the intensified uplift of the Tibetan-Himalayan complex, appears to have played an important role in Phrynocephalus diversification and evolution; (ii) intensification of aridity and geologic events in the Plio-Pleistocene and Quaternary glacial cycling probably had a significant influence on intraspecific diversification patterns within Phrynocephalus. Melville et al. [93] generated age estimates using two data. One method relied on mtDNA data, while the other relied on nuclear DNA data. The mtDNA data suggests that the estimated age of the common ancestor for Phrynocephalus was at 28.9 million years ago (Ma) with a 95% credibility interval of 36.2-21.1 Ma, a finding that is consistent with previous estimates [88]. Using the nuclear gene RAG-1excluding P. interscapularis, the estimate for the age of the common ancestor was somewhat younger, 15.8 Ma (95% CI 23–11.8 Ma), in accord with estimates produced by Guo and Wang [92]. However, both the analysis of Guo and Wang [92] and RAG-1 analysis did not include the basal species P. interscapularis, which was included in the mtDNA analysis by Melville et al. [93]. From this analysis, Melville et al. [93] concluded that including P. interscapularis, the age of Phrynocephalus would probably be Late Oligocene or Early Miocene, which coincides with the onset of aridification in Central Asia. Thus, their findings are consistent with a scenario that Phryncephalus is an old radiation, rather than it being a more recent diversification associated with the global deterioration of the Cenozoic climates when aridification intensified in the Late Miocene and into the Plio-Pleistocene. In addition, the deep phylogenetic split between the Chinese and Central Asian lineages of the P. alpherakii/P. guttatus/P. versicolor species group is probably related to the uplift of the Altay and Tien-Shan Mountains [93]. It was also suggested that gradual intensification of orogenetic processes and caused by them progressive aridization in the Middle East and Central Asia served as one of the factors favoring the evolution and initial diversification of the species complex of P. helioscopus [94].

Dunayev et al. [95] pursued the evolutionary history of 35 Prhynocepahlus species from COI segments (600-750 bp), which is, so far, the phylogenetic study with the most comprehensive taxonomic coverage. This phylogeny shows that the clade of Iranian species (e.g. P. scutellatus and P. maculatus) is sister to the remaining taxa, which joins Middle-Asian and Central-Asian taxa. Within the latter clade the representatives of the subgenus Megalochilus form a monophyletic group (P. mystaceus, P. interscapularis, P. sogdianus) of probably Irano-Turanian origin. The monophyletic Tibeto-Himalayan subgenus Oreosaura joins viviparous species (P. vlangalii, P. thebaldi, P. forsythia etc.) and is grouped together with oviparous species of Central Asia. Within Middle-Asian taxa three main species groups can be distinguished: (i) the groups of sun-watcher agamas (P. helioscopus, P. ersicus), (ii) P. ocellatus-P. strauchi group and (iii) P. guttatus-P. versicolor group. They further inferred that various bigoegoraphic events might cause diversification of Phrynocepahlus in continental Asia: starting with orogenetic processes on the territory of present-day Iranian and Tibetan plateaus, transgressions of Paratethys as major factors for old-splits, and up to aridisation processes and dynamics of major river valleys as factors determining recent speciation processes.

The topographic variation of the Qinghai-Tibetan Plateau, coupled with cyclical climatic changes in the Pleistocene, and alternating glacial-interglacial periods exerted the greatest influence on the current spatial distribution and genetic structure of viviparous toad-headed lizards. Jin et al. [96] studied the phylogeographic patterns of the Qinghai toad-headed lizard Phrynocephalus vlangalii by analysing sequence data from Mitochondrial DNA (mtDNA) sequences (partial ND2, tRNA^Trp and partial tRNA^Ala). The Qinghai toad-headed lizard is viviparous and lives at elevations of 2000–4500 m, and is the dominant terrestrial species in the QTP [26]. It is primarily distributed in the Qaidam Basin and surrounding area, with three major mountain chains, the Kunlun, the Arjin and the A’nyemaqen Mountains, in its range. With a thorough sampling across the entire distribution and a traditional phylogenetic analysis, five deeply diverged mtDNA lineages were recovered and geographical distribution of these lineages had almost no overlap [96]. Not only the spatial locations of the mountain chains coincided with the divisions of the lineages, but the times of the mountain uplift were also concordant with the times of mtDNA lineage divergence. Their NCPA analysis further demonstrated that the allopatric fragmentation resulted from historical vicariant events dominated the history of this species. Consequently, Jin et al. [96] concluded that the uplifts of the three major mountain chains formed the physical barrier that caused the initial vicariant events. They also inferred that populations from the Qaidam Basin appeared to have undergone major demographic and range expansions in the early Pleistocene (1.6 Ma), consistent with the colonization of areas previously covered by the huge Qaidam palaeolake, which desiccated at the onset of the Pleistocene [96].

Using data from 11 microsatellite DNA loci, Wang et al. [97] further supported the hypotheis of Jin et al. [96] that the uplift of the Arjin and the A’nyemaqen Mountains caused the initial vicariant events that led to the formation of the diverged lineages within P. vlangalii. Later, using a mitochondrial fragment ND4-tRNA^Leu, Guo et al. [98]confirmed that the uplift of the A’nyemaqen Mountains and glaciations since the mid-late Pleistocene, especially during the Kunlun Glaciation, are considered to have promoted the allopartric divergence of P. vlangalii. The diversification of P. putjatia may be triggered by the tectonic movement in the Huangshui River valley during the C phase of Qingzang Movement. Subsequently, the glacial climate throughout the Pleistocene may have continued to impede the gene flow of P. putjatia, eventually resulting in the genetic divergence of P. putjatia in the allopatric regions.

Jin and Liu [99] described the phylogeography of a unique endemic agamid lizard, Phrynocephalus erythrurus from the Qiangtang Plateau by analyzing sequence data from a mtDNA segment ND2- tRNA^Ala. P. erythrurus is viviparous and occupies the highest regions of any reptile on earth [26]. They found that (i) this species diverged into two major lineages/subspecies at 3.7 Ma corresponding to the Northern and Southern Qiangtang Plateau; (ii) the Northern Qiangtang lineage diverged into two subpopulations at 2.8 Ma separated by the Beilu River Region and Wulanwula Mountains [99]. Their NCPA analysis further demonstrated that the allopatric fragmentation and restricted gene flow were the most likely mechanisms of population differentiation. Their results also indicated the presence of at least three refugia since the Hongya glaciation. Consequently, Jin and Liu [99] concluded that the uplift of Tanggula Mountains movement and glaciations since mid-Pliocene have shaped phylogenetic patterns of P. erythrurus.

An aridification of the Tarim Basin and adjacent areas since middle Pleistocene has produced significant genetic structuring of the local fauna. Zhang et al. [100] compared the phylogeographic patterns, population structure and history of Phrynocephalus axillaris and P. forsythii using a mitochondrial fragment ND4-tRNA^Leu. They demonstrated that the two species might have experienced different evolutionary history throughout their current distribution. For P. forsythii, a vicariant event, as a consequence of geological isolation by the initiation of Quaternary folding at 2.1–1.2 Ma in the eastern Tien-Shan [101], and desert expansion, might have produced the significant divergence between the Tarim and the Yanqi populations. For P. axillaris, populations of the Yanqi, Turpan and Hami Basins might have been established through dispersal during demographic expansion. Climatic fluctuations caused alternate expansion and shrinkage of rivers and oases several times, which likely led to habitat fragmentation for both species. Interaction between vicariance, dispersal and habitat fragmentation produced the current distribution and genetic diversity.

4.2. Racerunner lizards (Eremias)

The lacertid genus Eremias Fitzinger 1834, is considered to comprise approximately 34 species, which inhabit sand, steppe, and desert regions from northern China, Mongolia, Korea, Central and Southwest Asia to Southeastern Europe. The reproductive biology of Eremias is notable in that there exist two reproductive modes: viviparity and oviparity. Most are oviparous, whereas the Eremias multiocellata complex (comprising 6 subspecies), E. buechneri, E. kokshaaliensis, E. yarkandensis, E. quadrifrons, and E. przewalskii are viviparous [102]. Despite the progress that has been made in recent years for the systematics of certain Eremias groups, the large-scale pattern of their evolution in time and space remains open. Here, we review the origin and diversification of racerunner lizards with their relationships to geological and climatic changes.

Until recently, our understanding of the phylogenetic relationships and historical biogeography of Eremias was mainly based on anatomical studies [103-105]. Boulenger [106] and FitzSimons [107] assigned most of the species now placed in Pedioplanis to the subgenus Mesalina within the large genus Eremias. Szczerbak [108] regarded Eremias sensu lato polyphyletic and considered Eremias (s. s.) endemic to Asia. Based on morphological characters and geographic distribution, Szczerbak [103] subdivided the inclusive genus Eremias (s. l.) into two distinct genera: the genus Mesalina as a north African and lowland Southwest Asian clade, and the genus Eremias (s. s.), which is endemic to Asia. Furthermore, Szczerbak [103] subdivided Eremias into five distinct subgenera: Eremias Fitzinger in Wiegmann, 1834 (group E. velox), Rhabderemias Lantz, 1928 (group Eremias scripta–Eremias lineolata), Ommateremias Lantz, 1928 (group E. arguta), Scapteira Fitzinger in Wiegmann, 1834 (group Eremias grammica), and Pareremias Szczerbak, 1973 (group E. multiocellata). The five subgenera were supported by Arnold (1986) on the basis of the hemipenial characters.

Wan et al. [109] made the first attempt to elucidate the phylogenetic relationships among the Chinese racerunner lizards on the basis of mitochondrial 16S rRNA data. They found that E. brenchleyi and E. argus formed a clade as the sister group of E. multiocellata. They further inferred that E. arguta, E.grammica and E. velox originated from Central Asia and the rest species in China originated from East Asia. However, these conclusions were tentative due to both limited taxa sampling as well as use of relatively short, partly nondiagnostic, gene fragment. On the other hand, one potential shortcoming of the analyses of Wan et al. [109] was that they did not incorporate secondary structural constraints of 16S rRNA into analyses.

Orlova et al. [110] studied the phylogeographic pattern of the steppe-runner lizard Eremias arguta by analyzing mtDNA cyt b gene sequences. The steppe-runner lizard inhabits steppes and semi-deserts of Eastern Europe and Middle Asia from Romania to Western Mongolia and China. They found an old split between E. arguta uzbekistanica and all other taxa, probably coming from the area of Ustyurt plateau. Consequently, Orlova et al. [110] inferred that this vicariant event is likely to have been caused by Paratethys regression.

More recently, a study [111 using combined DNA data sets (3925 bp) from two mitochondrial genes (cyt b, 12S rRNA) and one nuclear gene (RAG-1) and comprehensive taxonomic coverage resulted in three well-supported conclusions: (i) the species of the traditional genus Eremias form six clades; (ii) the Iranian plateau is the center of origin of the genus as a whole; (iii) genus Eremias should be divided into 4 distinct genera, and at least 14 new species should be described within these genera. Using a molecular clock calibrated with well dated paleological events as calibration, the average evolutionary rate for the mitochondrial genes was estimated as 1.3%-1.6% per million years. Based on the mitochondrial DNA data and using the globe molecular clock method, the origin time for the genus Eremias was dated to the early Miocene and the start of intrageneric divergence to the middle Miocene.

Another phylogenetic study [112] using mitochondrial 16S rRNA segment results in five major conclusions: (i) monophyly of Eremias and a clade comprising Eremias, Acanthodactylus and Latastia are recognized; (ii) monophyly of the subgenus Pareremias is corroborated, with Eremias argus being the sister taxon to Eremias brenchleyi; (iii) the first evidence that viviparous species form a monophyletic group was presented; (iv) Eremias probably diversified at about 9.9 million years ago (with the 95% credibility interval ranging from 7.6 to 12 Ma); (v) specifically, the divergence time of the subgenus Pareremias was dated to about 6.3 million years ago (with the 95% credibility interval ranging from 5.3 to 8.5 Ma), suggesting that the diversification of this subgenus might be correlated with the evolution of an East Asian monsoon climate triggered by the rapid uplift of the Tibetan Plateau approximately 8 Ma. A recent phylogeographic study [113] indicated that the Yellow River and Taihang Mountains may have acted as important barriers to gene flow in Eremias brenchleyi.

Rastegar-Pouyani et al. [114] explored the genealogical relationships and intraspecific differentiation of the rapid fringe-toed lizard Eremias velox by analyzing mtDNA cyt b and 12S rRNA gene sequences. The rapid fringe-toed lizard is widely distributed in the Iranian Plateau and Central Asia. In combination with molecular dating, they inferred that the E. velox complex originated on the Iranian Plateau in the Middle Miocene, which is in accordant with the result of Rastegar-Pouyani [115]. Another founding was that the northern Iranian clade diverged first some 11–10 Ma, caused by first uplifting of the Elburz Mountains in the late Miocene, and that the Central Asian lineages split from the northeastern Iranian lineage approximately 6 Ma, most likely as a result of uplifting of the Kopet-Dagh Mountains in the northern margin of the Iranian Plateau.

4.3. Conclusion

As exemplified with toad-headed agamas (genus Phrynocephalus) and racerunner lizards (genus Eremias), geomorphic and climatic changes in this big area definitely have remarkable influences on the regional and adjacent biogeographic patterns. Future comparative studies between Eremias, Phrynocephalus and other groups that distribute at the same area and inhabit similar habitat, can help to elucidate the origination, diversification and dispersion of these lizards, and test the link between speciation, adaptation, and historical changes in their biogeography. Specifically, further studies could search for differences in the ecological niches of species and clusters in order to assess ecological divergence among the several groups [116].

Appendix S1

List of the 422 squamates species.

Lacertilia:

Acanthosaura armata, A. lepidogaster, Alsophylax pipiens, A. przewalskii, Asymblepharus alaicus, A. ladacensis, A. sikimmensis, Ateuchosaurus chinensis, Calotes emma, C. medogensis, C. mystaceus, C. versicolor, Cyrtodactylus khasiensis, C. zhaoermii, Cyrtopodion dadunense, C. elongatum, C. medogense, C. stoliczkai, C. tibetanus, Dibamus bogadeki, D. bourreti, Dopasia gracilis, D. hainanensis, D. harti, D. ludovici, Draco blanfordii, D. maculatus, Eremias argus, E. arguta, E. brenchleyi, E. buechneri, E. grammica, E. kokshaaliensis, E. multiocellata, E. przewalskii, E. quadrifrons, E. stummeri, E. velox, E. vermiculata, E. yarkandensis, Eutropis longicaudata, E. multicarinata, E. multifasciata, Gehyra mutilate, Gekko auriverrucosus, G. chinensis, G. gecko, G. hokouensis, G. japonicas, G. kikuchii, G. melli, G. palmatus, G. reevesii, G. scabridus, G. similignum, G. swinhonis, G. taibaiensis, G. wenxianensis, Goniurosaurus bawanglingensis, G. hainanensis, G. lichtenfelderi, G. luii, G. yingdeensis, Hemidactylus aquilonius, H. bowringii, H. frenatus, H. garnotii, H. platyurus, H. stejnegeri, Hemiphyllodactylus typus, H. yunnanensis, Japalura andersoniana, J. batangensis, J. brevipes, J. brevicauda, J. dymondi, J. fasciata, J. flaviceps, J. graham, J. luei, J. kumaonensis, J. makii, J. micangshanensis, J. polygonata, J. splendida, J. swinhonis, J. tricarinata, J. varcoae, J. yulongensis, J. yunnanensis, J. zhaoermii, Lacerta agilis, Lamprolepis smaragdina, Laudakia badakhshana, L. himalayana, L. papenfussi, L. sacra, L. stoliczkana, L. tuberculata, L. wui, Leiolepis guttata, L. reevesii, Lepidodactylus lugubris, L. yami, Lygosoma bowringii, L. quadrupes, Mediodactylus russowii, Oriocalotes Paulus, Phrynocephalus axillaris, P. forsythia, P. guinanensis, P. guttatus, P. helioscopus, P. mystaceus, P. przewalskii, P. putjatai, P. roborowskii, P. theobaldi, P. versicolor, P. vlangalii, Physignathus cocincinus, Plestiodon capito, P. chinensis, P. elegans, P. liui, P. marginatus, P. popei, P. quadrilineatus, P. tamdaoensis, P. tunganus, Pseudocalotes brevipes, P. kakhienensis, P. kingdonwardi, P. microlepis, Ptyctolaemus gularis, Scincella barbouri, S. doriae, S. formosensis, S. huanrenensis, S. modesta, S. monticola, S. potanini, S. przewalskii, S. reevesii, S. schmidti, S. tsinlingensis, Shinisaurus crocodilurus, Sphenomorphus courcyanum, S. incognitus, S. indicus, S. maculatus, S. taiwanensis, Takydromus amurensis, T. formosanus, T. hsuehshanensis, T. intermedius, T. kuehnei, T. luyeanus, T. sauteri, T. septentrionalis, T. stejnegeri, T. sexlineatus, T. sylvaticus, T. viridipunctatus, T. wolteri, Teratoscincus przewalskii, T. roborowskii, T. scincus, T. toksunicus, Trapelus agilis, T. sanguinolentus, Tropidophorus berdmorei, T. guangxiensis, T. hainanus, T. sinicus, Varanus bengalensis, V. salvator, Zootoca vivipara.

Serpentes

Acalyptophis peronei, Achalinus ater, A. formosanus, A. hainanus, A. jinggangensis, A. meiguensis, A. niger, A. rufescens, A. spinalis, Acrochordus granulatus, Ahaetulla prasina, Aipysurus eydouxii, Amphiesma atemporale, A. bitaeniatum, A. boulengeri, A. craspedogaster, A. johannis, A. khasiense, A. metusia, A. miyajimae, A. modestum, A. octolineatum, A. optatum, A. parallelum, A. platyceps, A. popei, A. sauteri, A. stolatu, A. venningi, A. vibakari, Amphiesmoides ornaticeps, Archelaphe bella, Atretium yunnanensis, Azemiops feae, Blythia reticulata, Boiga cyanea, B. guangxiensis, B. kraepelini, B. multomaculata, B. nigriceps, Bungarus bungaroides, B. fasciatus, B. multicinctus, Calamaria buchi, C. pavimentata, C. septentrionalis, C. yunnanensis, Calliophis maculiceps, Chitulia inornata, C. ornate, C. torquata, Chrysopelea ornate, Coelognathus radiates, Cyclophiops doriae, C. major, C. multicinctus, C. ruffus, Daboia siamensis, Deinagkistrodon acutus, Dendrelaphis biloreatus, D.hollinrakei, D. ngansonensis, D. pictus, D. subocularis, Dinodon flavozonatum, D. rosozonatum, D. rufozonatum, D. septentrionalis, Disteira stokesii, Elaphe anomala, E. bimaculata, E. carinata, E. davidi, E. dione, E. schrenckii, E. zoigeensis, Emydocephalus ijimae, Enhydris bennettii, E. bocourti, E. chinensis, E.enhydris, E. plumbea, Eryx miliaris, E. tataricus, Euprepiophis mandarinus, E. perlacea, Gloydius blomhoffii, G. brevicaudus, G. halys, G. intermedius, G. lijianlii, G. monticola, G. saxatilis, G. shedaoensis, G. strauchi, G. ussuriensis, Hemorrhois ravergieri, Hierophis spinalis, Hydrophis fasciatus, H. gracilis, H. parviceps, Kerilia jerdonii, Lapemis curtus, L. hardwickii, L. colubrine, L. laticaudata, Leioselasma cyanocincta, L. melanocephala, L. spiralis, Liopeltis frenatus, Lycodon capucinus, L. fasciatus, L. futsingensis, L. gongshan, L. laoensis, L. liuchengchaoi, L. ruhstrati, L. subcinctus, L. synaptor, Macropisthodon rudis, Naja atra, N. kaouthia, Natrix natrix, N. tessellate, Oligodon albocinctus, O. catenatus, O. chinensis, O. cinereus, O. cyclurus, O. eberhardti, O. formosanus, O. lacroixi, O. lungshenensis, O. melanozonatus, O. multizonatus, O. ningshaanensis, O. ornatus, O. taeniatus, Oocatochus rufodorsatus, Ophiophagus Hannah, Opisthotropis andersonii, O. balteata, O. cheni, O. guangxiensis, O. jacobi, O. kuatunensis, O. lateralis, O. latouchii, O. maculosa, O. maxwelli, Oreocryptophis porphyraceus, Orthriophis hodgsoni, O. moellendorffi, O. taeniurus, Ovophis monticola, O. tonkinensis, O. zayuensis, Paratapinophis praemaxillaris, Pareas boulengeri, P. carinatus, P. chinensis, P. formosensis, P. hamptoni, P. margaritophorus, P. monticola, P. nigriceps, P. stanleyi, Pelamis platura, Plagiopholis blakewayi, P. nuchalis, P. styani, P. unipostocularis, Polyodontognathus caerulescens, Praescutata viperina, Protobothrops cornutus, P. jerdonii, P. kaulbacki, P. mangshanensis, P. maolanensis, P. mucrosquamatus, P. xiangchengensis, Psammodynastes pulverulentus, Psammophis lineolatus, Pseudolaticauda semifasciata, Pseudoxenodon bambusicola, P. karlschmidti, P. macrops, P. stejnegeri, Ptyas carinata, P. dhumnades, P. korros, P. mucosa, P. nigromarginata, Python bivittatus, Ramphotyphlops albiceps, R. braminus, R. lineatus, Rhabdophis adleri, R. callichroma, R. chrysargos, R. himalayanus, R. leonardi, R. nigrocinctus, R. nuchalis, R. subminiatus, R. swinhonis, R. tigrinus, Rhabdops bicolor, Rhadinophis frenatus, R. prasinus, Rhynchophis boulengeri, Sibynophis chinensis, S. collaris, Sinomicrurus hatori, S. kelloggi, S. macclellandi, Sinonatrix aequifasciata, S. annularis, S. percarinata, S. yunnanensis, Thermophis baileyi, T. zhaoermii, Trachischium monticola, T. tenuiceps, Trimeresurus albolabris, T. gracilis, T. gramineus, T. gumprechti, T. medoensis, T. sichuanensis, T. stejnegeri, T. tibetanus, T. yunnanensis, Typhlops diardii, T. koshunensis, T. lazelli, Vipera berus, V. renardi, V. sachalinensis, V. ursinii, Xenochrophis flavipunctatus, X. piscator, X. hainanensis, X. unicolor.

Acknowledgments

Support was provided by National Science Foundation of China (30700062, 31272281) and Knowledge Innovation Program of the Chinese Academy of Sciences (KSCX2-EW-Q-6).