Response of Distribution Range Against Climate Change and Habitat Preference of Four National Protected Diploderma Species in Tibetan Plateau

Lin SHI ,Xiudong SHI ,Yuning CAO ,Yayong WU ,Haijuan WEI ,Youhua CHEN ,Ziyan LIAO＊ and Yin QI,5＊

1 Chengdu Institute of Biology,Chinese Academy of Sciences,Chengdu 610041,Sichuan,China

2 University of Chinese Academy of Sciences,Beijing 101408,China

3 Yibin University,Yibin 644007,Sichuan,China

4 Tibet Autonomous Region Remote Sensing Monitoring Center for Ecological Environment,Lhasa 850000,Xizang,China

5 Mangkang Ecological Station,Tibet Ecological Safety Monitor Network,Changdu 854500,Xizang,China

Abstract Understanding the spatial distribution and habitat preference for rare and endangered species are essential for effective conservation practice.We examined the spatial distribution and habitat preference of four Diploderma species (Diploderma drukdaypo, D. laeviventre, D. batangense,and D. vela),which are endemic to the Qinghai-Tibet Plateau and are currently under state protection.We used the ensembles of small models (ESM) approach and predicted potential distribution ranges of the species in current and two future climate scenarios (SSP126 and SSP585).The degree of overlap between the predicted distribution ranges and existing natural reserves was further analyzed.Habitat preference was examined using a paired quadrat method.Our results predicted that D. drukdaypo has a current distribution range of 600 km2,which would decrease to 50 km2 and 55 km2 under the SSP126 and SSP585 respectively.For D.laeviventre,the current distribution range is 817 km2,with minimum changes in the two future climate scenarios (774 km2 and 902 km2).For D. batangense,the current distribution range is 875 km2,which would expand to 1 522 km2 and 3 340 km2 in the two future climate scenarios.Similarly for D.vela,the current distribution range is 1 369 km2,which would change to 1 825 km2 and 2 043 km2 respectively under the two future climate scenarios.The effect of protection of current nature reserves are likely low for those species;we found no overlap (D. drukdaypo, D. laeviventre) or little overlap(D. batangense 2.04%-3.56%, D. vela 15.52%-16.87%)between the currently designated protection area and distribution range under current and future climate scenarios.For habitat preference,stones appear to be the critical habitat element for those species although different species had different stone requirements.Taken together,we provided critical information on potential distribution ranges and habitat preference for four endangered Diploderma species,and confirmed the inadequacy of current nature reserves.The establishment of new or expansion of existing nature reserves is urgent for the conservation of those species.

Keywords climate change,distribution range,habitat preference,nature reserve,reptile

1.Introduction

With the increasing impact of global climate change,the suitable habitats for animals are no longer suitable,many species are experiencing distribution range reduction and habitats degradation (Ashrafzadehet al.,2019;Liet al.,2019;Wanget al.,2021).Understanding the spatial distribution and habitat preference for rare and endangered species is urgent for effective conservation practice(Morris,2003a;Rodríguezet al.,2007).

The distribution range reflects a long-term adaptation and likely experiences comprehensive influence of ecological factors (Ackerlyet al.,2010;Karamianiet al.,2018;Loiseauet al.,2020).With climate change,many species tend to migrate to high-altitude and high-latitude regions,resulting in changes in their distribution ranges(Midgleyet al.,2003).The size of distribution range can be used as an important indicator of threatened condition,and species with a wide distribution range may be less vulnerable to environmental changes(Sz?cset al.,2017).Analyzing the response of distribution range against climate change,and clarifying the spatial pattern of target species' distribution space will provide an important reference for effective protection on rare and endangered species(Satoet al.,2014).

Habitat preference significantly affects the distribution range,and protection of key habitat elements is the prerequisite for effective species protection (Ashrafiet al.,2013;Dela Torre and Nuneza,2021).Animals only can survive where their basic living requirements are met(Nemeset al.,2006).Habitat preference in animals is generally determined by environmental features that are favorable to animal survival,reproduction,and overall adaptability during long-term evolution and adaptation processes(Matthiopouloset al.,2015;Strauss,1991).It is an important process for animals to react and make decisions in response to available heterogeneous environments and has a certain inherent regularity.

Ecological niche models,also known as species distribution models(SDM),are commonly used to predict the actual and potential distribution of a species.The approach uses known distribution and related environmental data to construct models based on certain algorithms to determine the ecological requirements of species,and project the model to different times and spaces to predict the actual and potential distribution(Zhuet al.,2013).Modeling rare species is often difficult because they appear less frequently and many combinations of predictive variables can lead to overfitting of the model(Vaughan and Ormerod,2005).To overcome the overfitting issue in predictions on rare species,the ensemble of small model(ESM)has been proposed and successfully applied to predict small-ranged species(Breineret al.,2015).The ESM is a method that simulates species’ distribution by utilizing numerous small models,or simple bivariate models.These models consider only two predictors at a time and then combine all possible bivariate models to form an ensemble prediction.The advantage of this method greatly improves the ratio of the number of predicted variables to the number of observed values,which further improves the reliability of rare species distribution model predictions (Breineret al.,2018,2015),providing scientific reference for translocation protection,field investigation,and monitoring area selection for endemic or endangered species.

Dragon lizards of the genusDiplodermaare widely distributed in East Asia,particularly in southwestern China(Liuet al.,2023;Uetzet al.,2021;Wanget al.,2022).Of the 46 described species,32 are distributed in the Hengduan Mountains and adjacent areas.With the construction of hydropower and road in this region,the survival of manyDiplodermaspecies has been severely challenged,and a large number of original habitats have been severely destroyed(Shiet al.,2022).Four species ofD.drukdaypo,D.laeviventre,D.batangenseandD.velaof the region are currently listed as species of second class protection in the 2021 National Key Protected Wildlife List (Figure 1).Nevertheless,studies on conservation ecology of thoseDiplodermaspecies are scarce.We do not have any information about the spatial distribution and habitat selection of each species.We also do not know the conservation effect of current conservation systems on those species.TheD.laeviventreandD.velawere even assessed as data deficient in the China Red List of Biodiversity in 2021.This situation seriously affects the formation of conservation plans on those species.

In the present study,we examined the distribution ranges and habitat preferences ofD.drukdaypo,D.laeviventre,D.batangense,andD.vela.Our objectives are 1) predicting their distribution ranges under three scenarios,current and two future climate scenarios (SSP126 and SSP585),2) evaluating the degree of overlap between the current designated protected area and distribution ranges of each species,and 3)assessing habitat preference of each species using a paired quadrat method.This comprehensive investigation will provide valuable insights towards the conservation and management of these endangeredDiplodermaspecies,paving the way for informed decision-making and effective preservation efforts.

2.Materials and Methods

2.1.Study areaSpecies involved in this study are located in the dry and hot valley in the Hengduan Mountains of southeastern Qinghai-Tibet Plateau (see Figure 1).The fourDiplodermaspecies are distributed in the valleys of Nujiang,Lancang Jiang,and Jinsha Jiang,respectively,among whichD.laeviventreis distributed in the upper regions of Nujiang in Xizang,D.drukdaypois distributed in the valleys from Karuo to Chaya in upper regions of Lancang Jiang,D.velais distributed in the valleys from Mangkang-cuowa to Yanjing in the upper regions of Lancang Jiang in Xizang and also in the Lancang Jiang valley in Deqin County,northwestern Yunnan.D.batangenseis distributed in the Jinsha Jiang valley at the border of Mangkang County in Xizang and Batang County in Sichuan,as well as in the dry and hot valley of Batang County in the upper reaches of the Jinsha Jiang(Figure 1).

2.2.Prediction of distribution range with ecological niche modelsAs the first step of distribution prediction,we obtained the coordinates of actual distribution records of fourDiplodermaspecies through field survey (2020.06-2022.06),literature review and historical survey records (Floweret al.,2013;Thomaset al.,2012;Tsoaret al.,2007;Wanget al.,2019,2016).The specific information and sources are shown in Appendix 1.After removing duplicates and sampling points without determined geographic information,a total of 107 species occurrences were obtained,with 28 ofD.drukdaypo,26 ofD.laeviventre,14 ofD.batangense,and 39 ofD.vela.Overall,minimal number of occurrences meets the 10 X rule of thumb(Guisanet al.,2017)and can be used in ESM modeling(Liaoet al.,2022).

For model predictors,we used 19 bioclimatic variables,net primary productivity (NPP) and growing season length(GSL) from CHELSA V2 (CMIP6) dataset for both current(1981-2010) and future (2071-2100) (http://chelsa-climate.org/).The ecological meanings represented by different factors are shown in Appendix 2.We screened above 21 variables before modeling to avoid multicollinearity(Villemantet al.,2011) based on criteria with variance inflation factor (VIF) less than 5 using thevifstepfunction in the R package “usdm” (Naimi,2017).Finally,seven factors were selected including Bio2 (average of monthly diurnal temperature range),Bio4 (seasonality of temperature),Bio8(average temperature of wettest quarter),Bio13(precipitation of wettest month),Bio14(precipitation of driest month),Bio15(seasonality of precipitation)and GSL.

We predicted potential distribution range shifts under two climatic scenarios,a low-emission scenario SSP126 and a high-emission scenario SSP585 based on the IPCC Sixth Assessment Report (Kisslinget al.,2010;Meinshausenet al.,2020).For general circulation model,we choose the GFDL-ESM4 (Geophysical Fluid Dynamics Laboratory,Princeton,USA) model following the suggestion of the Intersectoral Impact Model Intercomparison Project(ISIMIP3b) (Dunneet al.,2020;O’Neillet al.,2016;Stoufferet al.,2017).

We used the ESM approach to predict the distribution range of theDiplodermaspecies under current and future climate conditions.We used four algorithms,namely Generalized Linear Model (GLM),Regression Tree (CTA),Artificial Neural Network (ANN) and Maximum Entropy Model (MaxEnt) because of their good ability to balance calculation time and prediction accuracy.The parameter settings of the four algorithms were adjusted by the bestfitting ESM,based on actual distribution data and 10 000 randomly obtained pseudo-absence points within the Qinghai-Tibet Plateau range(Breineret al.,2018).

We used the Somers'D score with 75%of the sampling data selected randomly for training and the remaining 25% used for testing.Bivariate models with a Somers'D score less than 0 were set to zero and were excluded from the construction of the ESM.The accuracy of the model predictions was evaluated based on area under receiver operating characteristic curve(AUC) and value of true skill statistic (TSS).The receiver operating characteristic curve(ROC)is a curve connected by all possible threshold points,with false positive rate as the abscissa and true positive rate as the ordinate (Phillipset al.,2006).The AUC value is the area of ROC curve enclosed by abscissas,which provides a single method of model performance assessment that does not depend on a specific selection threshold(Phillipset al.,2006).The AUC value ranges from 0 to 1,and the closer of the AUC value to 1,the better the predictive effect of the model.Models with as AUC value less than 0.7 likely have poor prediction accuracy,while models with an AUC value greater than 0.8 were proposed to have good prediction performance(Petersonet al.,2011;Zanket al.,2014).TSS is calculated by subtracting 1 from the sum of sensitivity and specificity,which tests the proportion of all validation points that are correctly predicted(Chuet al.,2017).The closer the TSS is to 1,the higher the prediction accuracy of the model.It is generally believed that the model prediction results can be used when the TSS is larger than 0.5.

We calculated the size of potential distribution area of each species based on binary predictions,and estimated their size changes under current and future climate conditions.To get the binary map,we used the maximum TSS threshold to convert the continuous habitat suitability into suitable(1)and not suitable (0).The potential area of current,the potential area of future,and the shrinking and extension area of future were calculated respectively.Importantly,based on several years of field sampling results and expert experience with the actual distribution of the four species (see section 2.1),we consider only the currently occupied ecological niches,specifically,we constrain the predicted distribution of the species to the currently found basins and their major branch streams.All calculations were carried out using ArcMap10.2(Esri,Redlands,CA,USA,http://www.esri.com).To evaluate the effectiveness of existing nature reserves in protectingDiploderma,we calculated the overlapping between the current and future distribution ranges of each species and existing natural reserves,and the proportion of protective area to the distribution range was calculated.

2.3.Habitat preference evaluationHabitat preference was measured using a paired quadrat method in the field following Shiet al.(2022).Specifically,the discovery point of a lizard was regarded as a selected site,and a 2 m×2 m quadrat was set around the selected site and used as a selected quadrat,while a control quadrat with similar size was determined by throwing a stone with random direction and distance.If a lizard was found in the random quadrat,we would discard it,and select a new control quadrat.The minimum distance between the random and selected quadrat was kept at 10 m to avoid spatial autocorrelation.We measured 11 factors from both selected and control quadrat,including habitat type,vegetation height,vegetation coverage,substrate size,stone size,stone coverage,distance to the road,slope aspect,slope,environmental temperature,and humidity.The definition on each habitat factor was listed in Appendix 3.Habitat type was determined by vegetation and topography,while vegetation height was measured with a tape measure at a precision of 1 cm.Vegetation coverage was visually estimated,and substrate size and stone size were measured using a tape measure at a precision of 1 cm.Stone coverage and distance from the road were visually estimated.Slope aspect and slope were measured using the Aovi app on a mobile phone.Environmental temperature was measured from a point of 10 cm above the ground with a handheld laser temperature gun (FR830,Shenzhen Jiakang Technology Co.,LTD),and humidity was measured using a handheld humidity meter(8813,Deli Group Co.LTD).

2.4.Statistical analysisThe habitat preference was analyzed by comparing habitat factors between selected and control quadrats(Wuet al.,2012).For categorical variables including habitat type,slope,and aspect,we used chi-square test.For vegetation coverage,rock coverage,and environmental humidity,we first transformed them into angular variables using inverse sine function,and then tested the distribution of all numeric habitat factors using S-W normality test.For habitat factors conformed to normal distribution,we used pairedt-test in comparing differences between selected quadrats and control quadrats.For habitat factors that did not follow normal distribution,the non-parametric Wilcoxon signed-rank test was used.All data analysis was conducted in R(version 3.6.3)(R Development Core Team,2019).

3.Results

3.1.Predictions of distribution rangeOverall,the ESM model performed very well in terms of prediction,with the minimum values of AUC and TSS for each algorithm being greater than 0.95(Appendix 4).ForD.drukdaypo,the potential distribution range is 600 km2under current climate conditions,mainly in the valley of Lancang Jiang from Changdu Karuo to Chaya section in Xizang.Under future climate change scenarios,the potential distribution range will greatly decrease.In the SSP126 scenario,the distribution range is approximately 50 km2,including 46 km2of current range,and 4 km2of new area,and a reduction of 554 km2.In the SSP585 scenario,the distribution is 55 km2,with 53 km2of current range,2 km2of new area,and a reduction of 547 km2.The distribution range ofD.drukdaypodoes not overlap with any nature reserves under current and future climate scenarios(Figure 2).

Figure 2 Distribution range of D.drukdaypo in current and future climate scenarios.A:range size of current climate conditions denoted by black colour;B: change of range size in future SSP126 climate scenario compared with current,with red polygon meaning reduced range,black polygon meaning invariant range,while purple polygon meaning increased range;C: change of range size in future SSP585 climate scenario compared with current.

ForD.laeviventre,the potential distribution range is 817 km2under current climate conditions,mainly in the valley of the upper region in Nu Jiang.The distribution range will not experience obvious change under future climate change scenarios.In the SSP126 scenario,the distribution range of this species is about 774 km2,with 649 km2of current range,125 km2of new area,and a reduction of 168 km2.In the SSP585 scenario,the distribution range will expand to about 902 km2,with 699 km2of current range,203 km2of new area,and a reduction of 118 km2.The distribution range of this species does not overlap with the nature reserves for current and future climate scenarios(Figure 3).

Figure 3 Distribution range of D.laeviventre in current and future climate scenarios.A:range size of current climate conditions denoted by black colour,while polygon of orange represents the potential nature reserve overlapped with range of D.laeviventre;B:change of range size in future SSP126 climate scenario compared with current,with red polygon meaning reduced range,black polygon meaning invariant range,while purple polygon meaning increased range;C: change of range size in future SSP585 climate scenario compared with current.

For theD.batangense,the potential distribution range is 875 km2under current climate conditions,mainly in the valley of Jinsha Jiang.The distribution range will experience obvious extension under future climate conditions.In the SSP126 scenario,the distribution range will be about 1 522 km2,including 807 km2of current area,715 km2of new area,and a reduction of 68 km2.In the SSP585 scenario,the distribution range is about 3 341 km2,including 2 724 km2of new area,mainly in the northwest,and a reduction of 258 km2.We found a little overlapping between distribution range and nature reserve inD.batangense,with 2.29% of area in Sichuan Zhubalong Provincial Nature Reserve in current climate conditions,in the future,theproportion of area in the nature reserve will be 2.04% and 3.56% in the SSP126 and SSP585 scenario,respectively(Figure 4).

Figure 4 Distribution range of D.batangense in current and future climate scenarios.A:range size of current climate conditions denoted by black colour,while polygon of orange represents the potential nature reserve overlapped with range of D.batangense;B:change of range size in future SSP126 climate scenario compared with current,with red polygon meaning reduced range,black polygon meaning invariant range,while purple polygon meaning increased range;C: change of range size in future SSP585 climate scenario compared with current.

ForD.vela,the potential distribution range is about 1 369 km2in current climate conditions.Under future climate conditions,the distribution range of this species will experience extension,mainly in the north.In the SSP126 scenario,the distribution range will be about 1 825 km2,with 1 223 km2of current area,602 km2of new area,and a reduction of 146 km2.In the SSP585 climate scenario,the distribution range will be about 2 042 km2,with 1 024 km2of current area,1 019 km2of new area in the north,and a reduction of 345 km2.We found a relatively high overlapping between distribution range and nature reserve inD.vela,with 16.87%of area in Mangkang Yunnan Golden Monkey National Nature Reserve current climate conditions,in the future,the proportion of protected area will be 15.95% and 15.52%in the SSP126 and SSP585 scenario,respectively(Figure 5).

Figure 5 Distribution range of D.vela in current and future climate scenarios.A:range size of current climate conditions denoted by black colour,while polygon of orange represents the potential nature reserve overlapped with range of D.vela;B:change of range size in future SSP126 climate scenario compared with current,with red polygon meaning reduced range,black polygon meaning invariant range,while purple polygon meaning increased range;C: change of range size in future SSP585 climate scenario compared with current.

3.2.Habitat preferenceForD.drukdaypo,we measured habitats from 50 quadrats,with 25 of selected quadrats and 25 of control quadrats.ForD.laeviventre,we measured habitats from 38 quadrats,with 19 of selected quadrats and 19 of control quadrats.ForD.batangense,we measured habitats from 44 quadrats,with 22 of selected quadrats and 22 of control quadrats.ForD.vela,we measured habitats from 96 quadrats,with 48 of selected quadrats and 48 of control quadrats.

Compared to control quadrats,D.drukdaypopreferred habitats with steeper slopes(X2=7.14,df=2,P=0.03),larger substrate sizes(control=35.56±5.32 cm,selection=83.6±13.77 cm,t=272,P=0.003),farther distances from roads(control=4.3±0.91 m,selection=4.92±0.96 m,t=230.5,P=0.020),and larger coverage of stones(control=58.96±4.93,selection=69.64±4.34,t=2.1374,P=0.043)(Figure 6).We found no significant difference in habitat type,slope direction,environmental temperature,environmental humidity,stone sizes,vegetation height,or vegetation coverage between selected and control quadrats(P＞0.05).

Figure 6 Habitat selection of four Diploderma species.A:D.drukdaypo;B:D.laeviventre;C:D.batangense;D:D.vela.SPS: substrate size;SC:stone covered area;DR:distance from road;SS:stone size;VH:vegetation height;Hu:environmental humidity.Gray represents the control group and black represents the select group.

D.laeviventrepreferred habitats with steeper slopes (X2=10.14,df=2,P=0.008),larger substrate sizes (control=43.74±7.54 cm,selection=136.89±16.98 cm,t=166,P＜0.001),further distance from roads(control=1.71±0.17 m,selection=3.97±0.55 m,t=4.72,P＜0.001),larger stone sizes(control=21.74±4.54 cm,selection=67.42±18.81 cm,t=152.5,P=0.022),and larger stone coverage area (control=41.84±4.73 cm,selection=76.89±7.66,t=172,P=0.002)(Figure 6).We found no significant difference in habitat type,slope direction,environmental temperature,environmental humidity,vegetation height,and vegetation cover between selected and control quadrats(P＞0.05).

D.batangensepreferred habitats with sunnier slopes (X2=6.02,df=1,P=0.014),higher environmental humidity(control=29.27±1.57,selection=36.82±3.19,t=52,P=0.009),larger stone sizes(control=25.08±9.52 cm,selection=46.8 ± 9.97 cm,t=4.0932,P=0.001),higher vegetation(control=67.69 ± 9.33 cm,selection=79.06 ± 9.83 cm,t=3.2529,P=0.009),and larger stone coverage area(control=32.11 ± 7.75 cm,selection=63.64 ± 5.91 cm,t=148,P=0.007) in scrub stone habitats (X2=7.09,df=2,P=0.024)(Figure 6).We found no significant difference in slope,environmental temperature,substrate sizes,vegetation cover,and distance from roads between selected and control quadrats(P＞0.05).

D.velapreferred habitats with larger stone sizes(control=18.49±1.9 cm,selection=28.2±3.28 cm,t=2.3068,P=0.026),higher vegetation (control=43.1 ± 13.09 cm,selection=57.88±6.77 cm,t=418,P=0.001),and larger stone coverage area(control=30±4.08 cm,selection=49.89±3.89,t=720,P＜0.001)in scrub stone habitats(X2=10.07,df=2,P=0.007)(Figure 6).We found no significant difference in slope,slope direction,environmental temperature,environmental humidity,substrate sizes,vegetation cover,and distance from roads(P＞0.05).

4.Discussion

In this study,we investigated the distribution ranges of four nationally protectedDiplodermaspecies in the Qinghai-Tibet Plateau under three climate condition as well as their habitat preferences.We also analyzed the overlapping between the distribution range of each species and existing nature reserves.

Climate change has an obvious impact on the distribution range ofDiplodermaspecies,but different species have different responses.The distribution range ofD.drukdaypois currently very small.Under two future climatic scenarios,its distribution range will greatly decrease.In contrast,the distribution range ofD.laeviventreis less affected by climate change,and its distribution range will not experience obvious change under future climate condition.The distribution ranges ofD.batangenseandD.velaare expected to increase under future climate conditions.The different responses may be resulted from different plasticity capacities toward climate change along elevation gradient (De Kortet al.,2016;Kingsolver and Buckley,2017).The range reduction inD.drukdaypois likely due to its loss of plasticity while adapting to high elevation through genetic assimilation,as tree speciesAlnum glutinosa(De Kortet al.,2016),and montane butterflyColias eriphyle(Kingsolver and Buckley,2017).Alternatively,the different responses ofDiplodermaspecies may result from their different geographical ranges,as suggested in CubanAnolislizards in which widespread species are more tolerant to changeable weather conditions,and can conserve their current geographical ranges under future climate change scenarios(Cajigaset al.,2023).

The distribution range shifts toward high altitude with climate change inDiplodermaspecies,although there is some variation between species.Under climate scenarios in the future,D.batangenseandD.velatend to expand their distribution ranges towards the northwest (Figures 4 and 5).The distribution range ofD.batangenseis expected to expand toward westward,whileD.velais expected to expand its distribution range toward several directions,with one tendency to the upper regions of Lancang Jiang.With global climate change,species distribution may shift to for favorable thermal environment,moving to higher altitudes is considered an important strategy to combat global warming (Sinervoet al.,2018) and has been found in many species,e.g.,Sin butterfly species in central Brazil (Wilsonet al.,2005) and amphibians in China(Duanet al.,2016).The previously unsuitable area would become suitable with global climate change in birds(Thomas and Lennon,1999) and the suitable area of Chinese giant salamander(Andrias davidianus)was predicted to shift towards high altitudes and latitudes in future climate conditions(Zhanget al.,2020).

The overlapping between the distribution range of eachDiplodermaspecies and current nature reserve is low.The distribution ranges ofD.drukdaypoandD.laeviventredo not overlap with the existing natural reserves in current and future climate conditions.There are some overlapping between the distribution range and nature reserves forD.batangenseandD.velain current and future climate conditions,but the overlapping ratio is low.ForD.batangense,the overlapping area in current climate conditions only makes up 2.29% of Sichuan Zhubalong Provincial Nature Reserve,in the future,the proportion of area in the nature reserve will be 2.04%and 3.56%in the SSP126 and SSP585 scenario,respectively.ForD.vela,the overlapping area in current climate conditions makes up 16.87% of Mangkang Yunnan Golden Monkey National Nature Reserve,in the future,the proportion of area in the nature reserve will be 15.95%and 15.52%in the SSP126 and SSP585 scenario,respectively.This is consistent with a recent meta-analysis on the conservation effect of amphibians and reptiles in Qinghai-Tibet Plateau by current conservation systems.There are more than 150 natural protected areas on Qinghai-Tibet Plateau,which covers 31.63%of the total area.However,over 70%of amphibian and reptile species are not receiving effective protection in the present,and over 300 amphibian and 500 reptile species will be completely lost due to climate change(Miet al.,2023).Therefore,conducting scientific investigations and evaluations,and optimizing the protection area system of the Qinghai-Tibet Plateau is very urgent.

Furthermore,we showed the habitat preference in fourDiplodermaspecies,which likely provide important information on future nature reserve design.D.drukdaypoandD.laeviventreprefers habitats with steep slopes,large substrate size,far from roads,and more rocks.D.batangenseselects habitats with sunny slopes,high vegetation height,more rocks and moderate humidity in shrub habitat,whileD.velaprefers habitats with shrub,high vegetation height and more rocks.On the one hand,the habitat selection may be related to differences in habitat composition,on the other hand,it may be related to species’ thermoregulation requirement (Morris,2003b).The habitat ofD.drukdaypois mostly composed of gravel,which likely helps them hiding in rock crevices and obtaining ample basking opportunities.Comparatively,D.laeviventremainly inhabits in bare rock zones with less vegetation,which likely allows them to receive more sunlight.Compared with gravel habitat inD.drukdaypo,D.laeviventreprefers a substrate with 140 cm of average rock radius,this may be related to large body size inD.laeviventre,and large rocks can provide shade for lizards during thermoregulation and anti-predation (Howes and Lougheed,2004;Monasterioet al.,2010).D.batangenseandD.velaprefer habitats with shrub,rocks,sunny slopes,which can facilitate basking and foraging simultaneously.Similar habitat preference was found inLacerta agilis,in which individuals preferred habitats with one meter of vegetation height(Nemeset al.,2006).The presence of shrubs or rocks can create a shading gradient for lizards to shuttle between sunlight and shade during body temperature regulation (Arnold,1987;Vachevaet al.,2020).More shrubs can provide abundant food resources for lizards because the availability of invertebrates in shrubs and grasslands is much higher than on rocks (Luppiet al.,2020).

In conclusion,our research provides important information on distribution range and habitat preference of four rare and endangered species ofDiploderma.We have also evaluated the sensitivity of differentDiplodermaspecies'distribution ranges against climate change,as well as the effectiveness of the existing natural protected areas forDiplodermaconservation.These findings provide direct evidence for delineation of protected areas and the critical protection elements for these four species ofDiploderma.Overall,the existing natural protected areas have not provided effective protection for the fourDiplodermaspecies,especially forD.drukdaypoandD.laeviventre.We recommend prioritizing the protection of theD.drukdaypo,followed byD.laeviventre,and timely delineating protected areas for these two species.A new nature reserve or key protection area should be established in the valley of Lancang Jiang from Changdu Karuo to Chaya section forD.drukdaypo,while the type locality ofD.laeviventrenearby Basu county should get the special protection.The conservation ofD.batangenseandD.velacan be improved through optimization of existing natural reserves by considering the northward migration under future climate change conditions.

AcknowledgementsThis work was supported by the Second Tibetan Plateau Scientific Expedition and Research Program (STEP) (Grant No.2019QZKK0402),Science and Technology Plan Project of Tibet Autonomous Region to Yin QI (Grant No.XZ202201ZY0030G),and China Biodiversity Observation Networks (Sino BON-Amphibian and Reptile).We are grateful to Qilong LIAO for his assistance in field investigations.Yin QI and Ziyan LIAO conceived and finalized the manuscript,Lin SHI and Xiudong SHI completed data analysis and prepared initial draft,Yayong WU contributed to the field survey.