Lizard richness in mainland China is more strongly correlated with energy and climatic stability than with diversification rates

Contemporary environmental, historical and evolutionary factors are increasingly used to decipher the drivers of spatial patterns of species richness. Evidence of such correlations for Chinese reptiles is scarce and poorly understood. We therefore explored the validity of the environmental capacity, historical climatic stability and diversification rates hypotheses on Chinese lizard richness.


| INTRODUC TI ON
Explaining patterns of species richness has long been a major focus of macroecological and biogeographic studies.At broad spatial and phylogenetic scales, species richness most often decreases with increasing latitude, a pattern known as 'the latitudinal diversity gradient' (Dobzhansky, 1950;Pianka, 1966).This pattern is prevalent in most major taxa, including birds (Belmaker & Jetz, 2015), mammals (Rolland, Condamine, Jiguet, & Morlon, 2014), reptiles (Roll et al., 2017) and plants (e.g.Qian et al., 2021).Some taxa, however, exhibit different patterns, including a bimodal gradient in New-World poeciliid fishes (see García-Andrade et al., 2021), and peak diversities away from the equator in Chinese insects, and New-World Lampropeltini snakes (Pyron & Burbrink, 2009).Notably, latitudinal gradients seem to vary more at smaller geographic (Liang et al., 2022) and lower level phylogenetic scales.The 'classical' gradient of diminishing richness with increasing distance to the equator becomes more prevalent in large clades studied at large geographic scales (Roll et al., 2017).
The environmental capacity hypothesis predicts that regions with higher ambient energy (i.e.temperature) (Hawkins et al., 2003), environmental productivity (Evans et al., 2005), and habitat heterogeneity (Hortal et al., 2013;Rahbek & Graves, 2001), could maintain more species (Barreto et al., 2021;Field et al., 2009;Hawkins et al., 2003;Hawkins & Pausas, 2004;Pianka, 1966;Storch et al., 2005).Correlations between richness and such factors, however, assume near instant responses (reflecting either speciation or, more likely, dispersal) to current-day conditions.Historical factors and evolutionary pathways (e.g.diversification rates, historical climates and habitat tracking), however, can greatly influence current diversity, through their effects on speciation, extinction and dispersal (e.g.García-Andrade et al., 2021;Miller & Román-Palacios, 2021).The historical climatic stability hypothesis posits that regions with stable climates will have more species because lineages have time to colonize, and speciate, in them (Haffer, 1969;Sandel et al., 2011).Stable climates may mean that inhabitants suffer lower extinction rates because taxa can adapt to the prevailing climate (Carolina & Moritz, 2008).Regions with less stable climates (across evolutionary time) are thought to have elevated extinction rates.The diversification-rate hypothesis envisions that speciation and extinction rates vary geographically (Jablonski et al., 2006;Tedesco et al., 2017).Regions (with similar ages and assuming equal dispersal) where net diversification rates are high (higher speciation or/ and lower extinction rates) will harbour more species (Miller & Román-Palacios, 2021;Rolland, Condamine, Jiguet, Morlon, & Moritz, 2014) if rates are maintained over long time periods.Species-specific diversification rates (also called tip diversification rate, DR) estimate present-day diversification rates, and are increasingly being used to explore spatial patterns of species richness (e.g.Jetz et al., 2012;Mitchell & Rabosky, 2017;Rabosky et al., 2018).That said, they could explain richness patterns only if such rates are consistent through both space and time, and dispersal rates do not vary.Drivers of richness, however, likely vary geographically (Velasco et al., 2018), temporally and between lineages (Pyron & Burbrink, 2009;Skeels et al., 2020).
Investigating the patterns and determinants for regional richness patterns across different taxa may, therefore, help to better understand broad processes generating local and regional diversity patterns.
While latitudinal diversity gradients are often depicted as showing smooth clines from the equator to the poles, deserts often harbour fewer species than expected given their latitudes (Fritz & Rahbek, 2012;Jenkins et al., 2013;Monteiro et al., 2022).Lizards, however, are well adapted to arid regions (Powney et al., 2010;Schall & Pianka, 1978), and exhibited richness patterns somewhat discordant to other vertebrates.Several studies (Lewin et al., 2016;Powney et al., 2010;Roll et al., 2017;Skeels et al., 2020) found that, while lizard diversity peaks in the tropics, similar to most other large clades, it is often also high in desert regions (though not as high as in the tropics).Thus, latitudinal richness gradients may actually be smoother in lizards than in many other clades.Studies of individual lineages could provide important insights for understanding how factors influence the distributions of those lineages.For example, each major lizard family in Africa has its own high-richness 'hotspots' which overlap little with hotspots of other taxa (Lewin et al., 2016).
The mechanisms driving richness patterns may likewise also vary within taxa and regions (Whiting & Fox, 2021).For example, temperature and topography account for richness patterns in most lizard groups, but net primary productivity is an important correlate of richness in some lineages (Skeels et al., 2020).
Reptile richness has been studied in regions which enjoy a wealth of research, such as North America (Schall & Pianka, 1978), Australia (Powney et al., 2010;Schall & Pianka, 1978), and Europe (Rodríguez et al., 2005), but is generally poorly studied elsewhere (but see Lewin et al., 2016;Tallowin et al., 2017 for Africa, New Guinea respectively), for example, in China.The reptile fauna of China is better known than that of many of its immediate neighbours (e.g.Laos, Myanmar, Kazakhstan and Kirgizstan) (http:// www.gardi nitia tive.org/).Continental China has highly variable climates, terrains and zoogeographical regions, over a wide range of latitudes (18.2°-53.6°),providing diverse ecological conditions for lizards to inhabit.At a slightly lower spatial scale (Figure 1), western China features mostly high elevations (the Qinghai-Tibet plateau, altitude: over 2600 m), while eastern China is mostly low (altitude: under 660 m).Middle longitudes are intermediate (mostly 660 ~ 2600 m), a pattern known as the 'three terrain stages' (Jiang & Yang, 2009, and Figure 1).Considering both environmental and climatic factors (e.g.temperature, precipitation, altitudes etc.),China is often divided into four major eco-geographical zones (but only three mainland ecoregions): the eastern monsoonal region, the northwestern dry region, the south-western Qinghai-Tibet plateau, and the southern tropics region (Liu & Shi, 2016).He et al. (2017) focused on the distribution terrestrial vertebrates at an even higher resolution.They proposed that China encompasses 10 major zoogeographic zones, and these are separated from each other based on the interplay of topography, geological context and current climate (Figure 1).Lizard clades often show strong F I G U R E 1 Sketch of analyses in this study.All analyses based on both grid cells (black solid lines) and terrestrial ecoregions (grey dashed lines, Olson et al., 2001) from the whole of China (all lizards and four major clades), the three terrain types (green), four ecogeographic regions (purple) and 10 zoogeographic regions (light blue).The three terrain types were based on the altitudes, where the western China is higher (over 2600 m) and the eastern China is lower (under 660 m).The four ecogeographic regions were based on the vegetation, climates, terrains, etc.The southern topics region is exclusive to islands -which are out of our study region (mainland China).While the zoogeography data are from He et al. (2017).T1-High stage; T2-Medium stage; T3-Low stage.E1-The south-western Qinghai-Tibet plateau (Tibet).E2-The northwestern dry region (Arid).E3-The eastern monsoonal region (East).Z1-East Himalaya.Z2-Tibetan Plateau.Z3-Northwest China.Z4-Northeast China.Z5-Inner Mongolia Plateau.Z6-Longzhong Plateau.Z7-North China.Z8-South China.Z9-Yungui Plateau.Z10-Taiwan, island regions were excluded from the analysis.elevation preferences.For example, the genus Gekko is mostly distributed below 2000 m (but see G. jinjiangensis, Hou et al., 2021) while agamids often have representatives at much higher altitudes (e.g.viviparous species like Phrynocephalus erythrurus can inhabit over 5000 m, Zhao et al., 1999).Lizards range across all Chinese latitudes and climates except the high alpine ones.At the same time, reptile biogeographic studies developed rapidly in China over the past 20 years which makes our knowledge of lizard distributions more comprehensive than in many neighbouring countries.For example, previous studies (Huang et al., 2011;Liang et al., 2022;Qian et al., 2007;Zhao et al., 2006) have explored patterns and ecological drivers of species richness among Chinese lizards and found strong environmental correlates of overall species richness.These studies, however, only focused on contemporary climates and did not examine richness relationships with evolutionary and historical factors.We aimed to fill this gap by generating new hypothesis of phylogenetic relationships for Chinese lizards.We aimed to fill this gap by generating new hypothesis of phylogenetic relationships for Chinese lizards, knowledge of which has recently been greatly revised and expanded (Liang et al., 2022;Wang et al., 2019).
We assessed whether current climatic and environmental factors (i.e.temperature, productivity, topographic heterogeneity), historical climatic stability and evolutionary factors explain lizard richness patterns in China.We test the following hypotheses: 1. Current environmental factors: We predicted that richness increases with mean annual temperature, and, to a lesser extent, with higher primary productivity and habitat heterogeneity, consistent with the environmental capacity hypothesis.
2. Historical climatic stability: We predicted the lizard richness would increase with increasingly stable climates.

Diversification rates:
We hypothesized that tip-based diversification rates would be positively correlated with richness.However, high tip diversification rates (observed over short time periods) may be associated with regions that have only recently been colonized by the lizard lineages that are rapidly diversifying there today.If so, it would result in negative (or no) relationships between diversification rates and richness.
4. We reasoned that the relationships between richness and the various predictors would be qualitatively similar across lizards in China overall and for each subclade.Likewise, we further predicted regional patterns to be similar to China-wide patterns, and the relationships at smaller spatial and phylogenetic scales to be generally weaker than those exhibited by all lizards across China.

| Species ranges
We collected comprehensive point-locality data for all 266 lizard species.We first compared these with published distribution data that include Chinese lizards (Liang et al., 2021;Roll et al., 2017) and amended them as necessary.We created new polygonal maps for species that had some point localities falling outside their published polygonal maps.For these maps, we followed the α-convex hull method, which requires at least three points and generates the lowest overestimation of range size (Meyer et al., 2017).We thus used new distribution maps for 125 lizard species, using a total of 1514 point localities (1-122 per species, mean: 12.95 ± 20.25; Appendix S2).Thirty-three species are known from only one or two locations in China.We buffered these locations with a 10-km diameter circle (area ≈ 78.5 km 2 ).Some species are distributed close to Chinese international borders.For these, we only kept the parts of the polygons that fall within China.For the other 141 species, we used maps from Roll et al. (2017) and Liang et al. (2021).The resulting geographical ranges for the 266 species of Chinese lizards can be found at Reptile Ecological Database of China (www.sinitic-repti le.top).
In order to avoid the potential effects of insularity on richness, we excluded all islands, and species endemic to them, from this study.Therefore, a total of 237 species were included in the following analyses.
A recent study revealed that drivers of Chinese lizard species richness patterns vary little with the scale of the grid cells examined (0.5°, 1°, 1.5°; Liang et al., 2022).The 100 km (approximately 1° × 1° at the equator) resolution of grid size is the one most commonly used in previous studies (e.g.Huang et al., 2011;Lewin et al., 2016;Rodríguez et al., 2005;Whiting & Fox, 2021).A finer resolution, while desirable, is impractical at the current state of our knowledge of the lizard fauna of some of the more poorly sampled regions of China.We excluded grid cells with less than 90% land area (i.e.partially at sea), or with less than 90% of their area within China (i.e. they are partially in bordering countries).We used the remaining 818 cells in our analyses.

| Phylogenetic tree
In order to generate most comprehensive phylogenetic tree for current Chinese lizards, we reconstructed a phylogenetic tree for 251 Chinese lizards by compiling DNA sequence data including seven mitochondrial (12S, 16S, CYTB, COI, ND1, ND2 and ND4), and two nuclear (BDNF, RAG1) genes, downloaded from the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov) and from a published thesis (Shi, 2010).
We aligned genes independently using MEGA 6 software (Tamura et al., 2013) using the default settings.We then concatenated all aligned sequences for the same species using 'phylotools' packages (Zhang, 2017).We performed phylogenetic reconstructions using both maximum likelihood (ML) and Bayesian inference (BI) approaches implemented in the PhyloSuite software (Zhang et al., 2020).We used ModelFinder to select the best-fit model using Bayesian information criterion (BIC) (Kalyaanamoorthy et al., 2017).
We imputed the phylogenetic position of the 15 species (see Appendix S3 for details) with geographic but no genetic data, using the SUNPLIN algorithm (Martins et al., 2013, https://bioin fo.inf.ufg.br/sunplin).We inserted the species missing from the phylogeny into the tree based on their taxonomy (i.e.placed them in their tree at a random place within their genera).We then generated 1000 such ML trees with all 266 species.Detailed information about GenBank accession number, trees can be found in Dryad (datad ryad.org,DOI: 10.5061/dryad.ht76hdrk3,also see Appendix S3 for a newick format phylogenetic tree).We dropped the island endemic species from the tree, and the following phylogenetic analyses (see below) were based on the 237 species inhabiting mainland China.

| Regional richness and richness of major clade scales
We obtained the spatial ranges of the 'three terrain stages', the 'four major eco-geographical zones', and '10 zoogeographic zones' (but excluding insular regions), from Qi et al. (2017), Liu and Shi (2016) and He et al. (2017), respectively.We explored patterns and drivers of richness within the 'three terrain stages' (i.e.high, middle, low elevations), within the four eco-geographical zones and within the mainland (without any islands) zoogeographic zones of He et al. (2017).
The unit of analysis for all data sets (China, three terrains, four ecogeographical zones and zoogeographic zones) was the 100 × 100 km grid cell.For mainland China overall (only) we also tested richness within the 63 terrestrial ecoregions of Olson et al. (2001) as the unit of analysis (see Figure 1), controlling for (log 10 -transformed) ecoregion area (see below).
We also tested our hypotheses within each of the four major lizard clades: Agamidae, Lacertidae, Scincidae and Gekkota (including members of the Eublepharidae, Gekkonidae and Sphaerodactylidae) that together comprise 229 species (97% of mainland lizard species).
Anguids, Dibamus, Shinisaurus, and varanids (1-2 species in each) were used in these analyses.The lineage-based analyses were only conducted for China overall with gridded richness as the response variable.Jetz et al. (2012) proposed a tip-based measurement of diversification rates (called tip rates) for extant species, that is, species-specific diversification rates (the inverse equal splits rate, see equation 1 for tip DR calculation).Such tip rates are considered to provide good representation of speciation rates (Miller & Román-Palacios, 2021;Title & Rabosky, 2019), because they capture only the recent purebirth diversification rate (Jetz & Fine, 2012;Title & Rabosky, 2019).

| Diversification rate
We estimated mean tip diversification rates (mDRs) among cooccurring species in each grid cell and terrestrial ecoregion.We then evaluated the relationships between species richness and mDRs to assess the potential effects of diversification rates on lizard richness.
Equation 1: calculation of mean tip diversification rates.N i is the number of edges on the path from species i to the root, and EL i is the length of internal branch j from the branch to the root (Cai et al., 2018;Jetz et al., 2012).
For the analysis of China's 63 terrestrial ecoregions, we included log 10transformed ecoregion area as an additional predictor to correct for species-area relationships (Thierry et al., 2011, also  to affect the environment in a multiplicative rather than a linear fashion (e.g. 100 mm more rain per year is negligible in tropical rainforests, but makes the difference between severe drought and a normal year in deserts, while a 15% increased precipitation likely have more similar effects in both regions).Transformation further helped to normalize model residuals.A negative correlation between climatic variability and richness would be consistent with higher climate stability driving higher richness.We extracted the values of these climatic and environmental factors for each grid cell and each of the 63 terrestrial ecoregions using the 'raster' package (Hijmans, 2020) in R (R core team, 2021).

| Statistical analyses
We analysed richness in grid cells (at all scales) and ecoregions (only for overall China scale) (Figure 1).We tested how species richness varies with latitude using generalized additive models (GAMs) for lizards overall in grid cells and the 63 terrestrial ecoregions.We tested the four major lizard clades separately using grid cell richness.GAM was chosen because it can capture non-linear relationships between predictors and the response.We did not use GAMs to test latitudinal gradients at lower scales (terrain, ecogeography, zoogeography), because each subregion has relatively narrow latitudinal range.
In order to test how contemporary environmental factors (temperature, NPP, altitude heterogeneity and area), historical factors (short-and long-term temperature and precipitation stability) and evolutionary factors (diversification rates) are related to richness, we used spatial autoregressive (SAR) models (Dormann et al., 2007), to account for spatial autocorrelation.We also calculated the variance inflation factor (VIF) for these predictors of each model to see if there is a multicollinearity.
We performed SAR models to assess if there is a difference between long-term and short-term past climates stability (see Table 1 for model details) in shaping current lizard richness.These We applied piecewise structural equation models (pSEM, Lefcheck, 2016) to jointly evaluate the proposed hypotheses, considering direct and indirect (Barreto et al., 2021) effects (Figure 2).
PSEM analyses were performed across China for lizards overall (all grid cells and 63 terrestrial ecoregions) and within the four major clades.Separate pSEM analyses using grid cell richness as the unit of analyses were run for the three topographical, four ecogeographical and 10 zoogeographical regions (for all lizards; Figure 1).In the sub-China analyses, we only analysed subregions with >15 grid cells (i.e. without islands in ecogeography, Taiwan and Longzhong Plateau for zoogeography, see Table 2).For pSEM analyses, we included LGM but not CSI climate stability in view of our SAR results (see results).
In order to assess the relative importance of each factor, we performed variance decomposition analyses [as described by Gross et al., 2017].These were done only for regressions exploring the direct relationship between species richness and predictors.

| Sensitivity analyses
Repeated species co-occurrences across multiple grids/regions can generate strong and spurious relationships between richness and the factors used to assess it (Hawkins et al., 2017).In order to detect such influences, we generated 100 random databases by randomizing raw lizard richness and mDR separately across grids and across the 63 terrestrial ecoregions.Then, we repeated SAR and GAM analyses and calculated the 100 randomized Nagelkerke pseudo-R 2 (for SAR) values (Hawkins et al., 2017).Only when all the 100 randomized values were significantly lower than the observed values did we consider the observed relationships reliable.

| Latitudinal and altitudinal patterns
We found that lizard richness was highest in southern China at with increasing latitude between 20° and 30° and more slowly further north (Figure 3b).Richness likewise decreased with increasing latitude in skinks, geckos and agamids, while lacertids showed no latitudinal gradients (Figure 4).

| Correlates of lizard richness
We found low collinearities between long-and short-term past climate stability indices (all VIF values less than 2).The SAR models (see Table 1) showed that lizard richness was more strongly correlated with the short-term past climate stability than with longterm stability.Lizard richness was negatively correlated with the variation in temperature (in most cases at both ecoregions and grid cells) and precipitation (in grid cells) (i.e.regions with more stable temperature or precipitation have more lizards) only since the LGM (Table 1).Richness was independent of longer term (CSI) stability in precipitation and temperatures in most cases at both the grid-cell and ecoregion scales (Table 1).
The structural equation models (pSEMs) showed that lizard richness overall, and within three of the four major clades, was consistently positively correlated with temperature (Figures 5 and 6).Elevational heterogeneity was positively correlated with lizard richness across the 63 ecoregions (Figure 5b) and in the Scincidae (Figure 6), and net primary productivity was positively correlated with species richness only in the Scincidae (Figure 6).Diversification rates were negatively correlated with lizard richness overall and within most major clades (Figures 5 and 6).Only Gekkota richness showed positive correlations with diversification rates (Figure 6).Temperature was also negatively correlated with diversification rates (i.e.indirect path-7 in Figure 2) across grid cells and in Agamidae and Lacertidae (Figures 5 and 6).
Lizard richness was negatively correlated with past short-term climates stability (since LGM) across the 63 ecoregions and the grid cells (Figure 5), and in the Scincidae (Figure 6).
At regional scales (Table 2), we found similar relationships to those found for lizards overall, and within major clades.This was especially true for temperature, which was positively correlated with lizard richness in most regional analysis (all terranes and ecogeographic regions, and 5 of 8 zoogeographic regions) (Table 2).Net primary productivity and elevational heterogeneity were positively correlated with lizard richness only in western (i.e.Tibet) and eastern China (low altitude regions) (Table 2).Diversification rates were negatively correlated with lizard richness in most regions.There were positive relationships between lizard richness and diversification rates only in three regions (East Himalaya, Inner Mongolia Plateau and North China zoogeographic regions) (Table 2).Finally, lizard richness was positively correlated with stable climates (temperature, precipitation or both temperature and precipitation) since the LGM in western and northern China (Table 2).
We found temperature has the highest importance in 11 of 20 models (Appendix S4), mDR explains the highest variances in 5 of 20 models but these correlations were negative in two of these five models.This is contrary to the diversification rates hypothesis that predicts a positive relationship (Appendix S4, Figures 5 and 6, Table 2).Hence, we do not consider the support for this hypothesis to be strong.

TA B L E 2
Piecewise structural equation modelling for drivers at regional scales.

| Sensitivity analyses
Observed Nagelkerke pseudo-R 2 values were larger than all random values for both species richness and mean diversification rates for lizards overall and within the four major clades, and at both the entire China and regional scales (see Appendix S5 for detailed tested factors and results).This suggests that our results are not merely statistical biases.

| DISCUSS ION
Understanding the mechanisms that generate and maintain patterns of species richness can provide valuable insights for understanding diversity (Cerezer, Cáceres, & Baselga, 2022;Cerezer, Machac, et al., 2022) and predicting the impacts of future climate changes on it.We described the geographical pattern of Chinese lizard richness in detail at the country and regional scales, and We did not include snakes in this study, because of worries regarding the less complete knowledge of their phylogeny and distribution (compared to lizards).Previous global studies of squamates (Caetano et al., 2022;Roll et al., 2017) found that snake richness follows a very strong latitudinal pattern in China, as well as a strong East (rich)-West (poor) gradient between latitudes ~29 and 36° N, with snakes being completely absent from most of the Tibetan Plateau.This agrees quite well with what is known for lizards as a whole (Figure 3a) but differs in some important details from patterns shown by major lizard (and probably snake) clades (Figure 4).Thus, we hypothesize that our general results, and hence our major conclusions, are robust to the inclusion or exclusion of snakes.

| The latitudinal patterns of lizard and major clade diversity
The latitudinal diversity gradient is a common pattern at both global (e.g.Hillebrand, 2004) and regional (e.g.reptiles in North America, Whiting & Fox, 2021) scales (but see García-Andrade et al., 2021 for fishes).Chinese lizards, as previously reported (Huang et al., 2011;Qian et al., 2007;Zhao et al., 2006), show richness patterns consistent with the latitudinal diversity gradient at both the grid and ecoregion scales (Figure 3).This may be because most lineages showed a latitudinal diversity gradient, although details differed (see Figure 4).thus weakening the effect.Likewise, the most topographically varied regions probably contain montane climates that are detrimental to most lizard species, lowering overall richness.Primary productivity may enable higher population densities potentially buffering against extinction (due to higher amounts of available energy, Figures 2,5) and allow the existence of more specialist species (Abrams, 1995;Srivastava & Lawton, 1998).However, they may also facilitate increased numbers of competitors, predators and pathogens, which can lower richness (Braz et al., 2020).
We found that lizards are very sensitive to current climate (especially to temperature; as was found many times in the past).We further show they are sensitive to climate change as expected by their ectothermic physiology.World climates have repeatedly changed, and the Pleistocene was an especially dynamic period, characterized by frequent and drastic temperature fluctuations.These likely strongly affected the distributions of lizards.It is therefore unsurprising that we found effects of mostly short-term (i.e.since LGM) climatic stability (Table 1).We suspect that climate changes since the LGM caused major shifts in lizard distributions (Herrando-Moraira et al., 2022).Such changes likely swamped, and might have erased, much of the legacy of older climatic changes.The findings that shortterm climate changes strongly affect lizard distributions suggest that the impacts of current climate changes (Araújo et al., 2008) on Chinese lizard biodiversity is likely to be severe.
We found negative relationships between richness and diversification rates in many cases.This counterintuitive result has been previously found in several taxa, including terrestrial vertebrates (Cerezer, Machac, et al., 2022;Hutter et al., 2017;Tejero-Cicuéndez et al., 2022) and fishes (Egan et al., 2022;Hortal et al., 2013).Such a negative relationship can stem from faster speciation being selected for in extreme environments, which do not allow enough time for biodiversity to accumulate (Harvey et al., 2020).Indeed, much time probably needs to pass for diversification rates to manifest in richness patterns (Miller & Román-Palacios, 2021), especially if dispersal in response to short-term climatic fluctuations is rife.An additional option is that the mDR metric fails to capture clade-wide species accumulation.But past climates worldwide, and in China specifically, have changed frequently (Herrando-Moraira et al., 2022).Therefore, we would not expect a strong positive correlation between diversification rates (that operate over longer time scales) and current species richness.Some additional explanations may apply to the case of Chinese lizards (and potentially other systems).A non-mutually exclusive explanation can be that such extreme environments see higher emigration rates, where newly evolving species disperse elsewhere.
We suggest that current Chinese lizard richness patterns (especially in western regions, e.g., Tibet) could be the outcome of dispersal to recolonize northern regions since the end of the last Glacial period.
South-eastern China provided a glacial shelter for many species from extinction, it was never glaciated and its climate was less influenced by Quaternary glaciations than other Chinese regions (Wang, 2008;Zhang, 2002).In summary, relative to the long periods needed for diversification (Miller & Román-Palacios, 2021), the short temporal scale of dispersal processes, due to recent climates change (relative to speciation times), may swamp the evolutionary signal.Furthermore, they could help explain why we found stronger relationships between species richness and current climatic and environmental factors than with historical ones.
The relationships between richness and the predictors across geographic scales, and within major clades.Elsewhere-for example in Africa, different diversification centres for different clades result in clade specific hotspots (Lewin et al., 2016).China, however, has far fewer species and lineages (~270 species from 10 families) than Africa.Furthermore, parts of China (the Tibetan Plateau and the North) are much colder than anywhere in Africa-and these are prohibitive for all lizard lineages.

| CON CLUS ION
Taken together, our results suggest that lizard richness patterns in mainland China followed current environmental conditions (environmental capacity hypothesis) and short-term historical stability (historical climatic stability hypothesis).These factors explain most contemporary variation in richness.The diversification rate hypothesis and long-term climatic stability hypothesis were unsupported.
We suspect that contemporary species richness patterns are often more strongly influenced by dispersal and local extinctions than by in situ diversification.
see results).In order to test whether climatic stability (at different timescales) affects lizard richness patterns, we used two different time periods (short and long term) of 21,000 and 3.3 million years, respectively, to explore the roles of historical climate stability (in both temperature and precipitation) in shaping lizard richness pattern.We posit that the longer timescale is potentially more influenced by speciation while the shorter one probably mostly reflect range shifts (dispersal) and extinctions.We used newly published historical climate stability data(Herrando-Moraira et al., 2022) as the proxy for long-term effects.These are the most comprehensive data available for the period ranging from the Pliocene (~3.3 Ma) to the present (Herrando-Moraira et al., 2022).Areas with low climatic fluctuations had low climate stability index (CSI) values.We also calculated climatic stability since the last (1) Equation 1: tip DR i = ∑ LGM; ∼21 ka) to today.This period includes several main climatic regimes including the LGM glaciations, warmer climates towards the end of the Pleistocene, the Younger Dryas Stadial (YDS, ~14 ka) episode of cooling and warmer Holocene climates.We downloaded three climate layers: 21 ka (LGM), 14 ka (YDS) and 1980-2010 (current), for the LGM, the YDS and the current climate, respectively, from climatologies at high resolution for the Earth's land surface areas (CHELSA, https://www.chelsa-clima te.org/.Karger et al., 2021), and rastered the mean values of each layer for each grid and ecoregion.We then calculated the standard deviation (SD) values of each mean climatic variable across the three periods for each grid cell and ecoregion (lower SD means higher climatic stability).Before calculating the SD, we log −10 transformed precipitation values since precipitation is likely models included (1) all factors, (2) all factors except long-term climates stability (i.e.only with LGM climate stability), (3) all factors except short-term climates stability (i.e.only with CSI climate stability), (4) all factors without past precipitation stability (i.e.only has past temperature stability), (5) all factors without past temperature stability (i.e.only has past precipitation stability).TA B L E 1 Results of multivariable spatial autoregressive (SAR) models.Analyses of all lizards for mainland China, and each column represents a SAR model of grid and ecoregion.All these six models are different combinations of long-and short-term historical climatic stability at grid and ecoregion scales, to see which factor plays a stronger role.Values are regression slope of each model, ***p < 0.001, **p < 0.01, *p < 0.1.NA means factor was not included in the model (column).AIC was the Akaike information criterion value of each model (column).Boldface values mean results were significantly different from the null at alpha = 0.05.Nagelkerke's pseudo-r 2 value was provided for each model, but these cannot be interpreted as the percentage of variance explained by the model.
both the grid and ecoregion scales, and richness was lowest in northeastern China and in the Qinghai-Tibet plateau in the west (Figure 3a,b).Lizard richness, overall, decreased sharply (Figure 3b) F I G U R E 2 Theoretical model structure and hypothesis tested in our spatially explicit structural equation models.The thick lines (1-3) represent our main hypotheses (direct relationships); line 4 represents the direct relationship of area via the species-area relationships.Lines 5-9 represent indirect relationships (see methods for details).Red lines: positive correlations, blue lines: negative correlations.SR = species richness, AltSD = altitude standard deviation, Tem = temperature, NPP = net primary productivity, mDR = mean diversification rates, HTSD = historical temperature standard deviation, HPSD = historical precipitation standard deviation since the LGM.

F
Latitudinal lizard richness patterns at the grid (a, b) and ecoregion (c, d) scales.Maps in equal-area Behrmann projection.b & d: mean and 95% confidence intervals (in grey) for richness at different latitudes based on generalized additive models.F I G U R E 4 Species richness (a, c, e, g) and latitudinal gradients (b, d, f, h) within the four major clades of Chinese lizards.Maps in equal-area Behrmann projection.Models of richness along latitude (b, d, f, h), and their 95% CI (grey) are based on generalized additive models.SR: species richness.Representative species within each group are; Phrynocephalus mystaceus (a): Agamidae; Teratoscincus scincus (c): Gekkota; Eremias arguta (e): Lacertidae; and Plestiodon chinensis (g): Scincidae (Photos: Tao Liang).explored, for the first time (to our knowledge), potential historical and evolutionary drivers.Grid-cell and ecoregion scales revealed different correlates of lizard richness.Ecoregions probably reflect dispersal limitations of species assemblages to some extent (i.e. they represent 'real' assemblages with meaningful ecological associations).Grids are completely arbitrary units in both size and boundary placement and may combine allopatric assemblages together.But analyses based on grids may present higher variance in the predictors and in richness.Furthermore, the grid cells we use are smaller than the size of even the smallest ecoregion.Hence, ecoregions are more likely to aggregate allopatric species-especially in regions with multiple small range endemics.This may have caused us to fail to detect the predicted correlations based on ecological and evolutionary frameworksdespite the higher sample sizes (i.e. higher statistical power) of the grid-based analyses.Current levels of ambient energy emerged as the strongest correlate of Chinese lizard richness patterns, while evolutionary and historical factors seemed less important.
Within the four major clades, lacertids differed from the other three lineages and do not exhibit a clear latitudinal pattern of species richness (Figure4).The two major lacertid genera in China, Eremias and Takydromus (each with 15 species, 47% of the Chinese Lacertidae), have the opposite latitude trends.This may reflect the results of historical dispersal processes, since this family originated and diversified at relatively high latitudes in the Western Palearctic and spread to eastern Asia later (Garcia-Porta et al., 2019).Indeed, the Palaearctic Eremias species are mainly found in arid northern China, and their diversity decreases southwards.The F I G U R E 5 Piecewise structural equation modelling for drivers of overall Chinese lizard richness (a: grid cells, b: 63 ecoregions).The pathways show how the factors are linked and drive the lizard richness.Values represent standardized coefficients (the thick of each line is correlated with the absolute value of standardized coefficient).Blue arrows denote significant negative correlations and red arrows denote significant positive correlations.SR = species richness, AltSD = altitude standard deviation, Tem = temperature, NPP = net primary productivity, mDR = mean diversification rates, LGM-T = historical temperature standard deviation since LGM, LGM-P = historical precipitation standard deviation since LGM.
Note: ***p < 0.001, **p < 0.01, *p < 0.1.Values are the standardized coefficients of each model.Historical climate stability including both temperature (left) and precipitation (right) stability.Boldface values mean results were significantly different from the null at alpha = 0.05.Heterogeneity is the standard deviation (SD) of altitudes of each grid.The details on regions are found in Figure1.