Variable vulnerability to climate change in New Zealand lizards

The primary drivers of species and population extirpations have been habitat loss, overexploitation and invasive species, but human‐mediated climate change is expected to be a major driver in future. To minimise biodiversity loss, conservation managers should identify species vulnerable to climate change and prioritise their protection. Here, we estimate climatic suitability for two species‐rich taxonomic groups, then use phylogenetic analyses to assess vulnerability to climate change.


| INTRODUC TI ON
Anthropogenically driven biodiversity loss is predicted to increase in the future (Pecl et al., 2017). The main drivers of species extinctions and population extirpations in the recent past have been habitat loss and degradation, overexploitation and introduction of invasive alien species (Dirzo et al., 2014). However, several lines of evidence suggest that human-mediated climate change will increasingly be a prominent, if not the leading, cause of extirpations in the coming century (Dirzo et al., 2014). To minimise future biodiversity loss, conservationists need to identify species likely to be vulnerable to the impacts of climate change (Pecl et al., 2017).
High-profile meta-analyses and reviews of the response of species to recent and future climate change have primarily focused on the Northern Hemisphere and the tropics (e.g., Chen et al., 2011;Parmesan & Yohe, 2003). Ectothermic vertebrates in these regions appear to be particularly susceptible to climate change (Huey et al., 2009(Huey et al., , 2012, with tropical vertebrates typically being predicted to be more vulnerable to increasing temperatures than their temperate counterparts (Huey et al., 2009(Huey et al., , 2012; but see, e.g., Sunday et al., 2014). The majority of climate change research has focused on ectotherms, such as lizards in the tropics (Huey et al., 2009(Huey et al., , 2012Sinervo et al., 2010). However, comparatively little is known about how climate change is likely to impact lizard distributions in temperate regions, and the Southern Hemisphere is particularly understudied.
In Aotearoa New Zealand (NZ), a temperate archipelago, current predictions of extinction risks from climate change are estimated to be high compared to other geographical regions (Urban, 2015). This is because NZ harbours diverse assemblages of endemic species with small ranges, and there are relatively small land areas that can limit shifts to new habitat (Urban, 2015). However, the threat of climate change has not been rigorously modelled for most NZ species, including lizards. The higher-level taxonomic diversity in NZ lizards is limited to diplodactylid geckos and eugongyline skinks (Chapple, 2016;Chapple et al., 2009;Nielsen et al., 2011), with a single highly speciose endemic radiation in each group. Including undescribed taxa there are at least 61 extant skink species and 43 extant gecko species, all of which are endemic (Hitchmough, Barr, et al., 2016;Hitchmough, Patterson, et al., 2016). This lizard diversity is remarkable for a temperate region (Chapple, 2016;Hitchmough, Barr, et al., 2016;Hitchmough, Patterson, et al., 2016), and many species have life-history traits such as large body size, low reproductive output and late maturity that may make them particularly vulnerable to climate change (Chapple et al., 2021;Hoffmann & Sgrò, 2011;Nelson et al., 2014;Sinervo et al., 2010;Whitaker, 1978). Numerous species have experienced substantial range reductions following human settlement and the introduction of mammalian predators, while others have apparently naturally restricted distributions (Hitchmough, Barr, et al., 2016;Hitchmough, Patterson, et al., 2016;Nelson et al., 2014). A recent study suggested that NZ lizards with a larger body size had a higher threat status, as did species with increased habitat specialisation and decreased range size . Larger bodied lizard species are well known to be more vulnerable to predation by invasive mammalian predators than smaller bodied species in NZ (because the refugia they use when they are inactive and unable to flee predators are larger and therefore more accessible to predators), resulting in disproportionately large range contractions in these species and hence potentially greater vulnerability to climate change (Hoare et al., 2007;Nelson et al., 2014;Whitaker, 1978). Thus, conservation efforts have been primarily aimed at control and eradication of invasive mammalian predators, as well as reintroductions to restore populations (Nelson et al., 2014;Towns et al., 2016).
In this study, we examine the vulnerability of NZ lizards to projected climate change. Predictive models such as correlative species distribution models (SDMs, often termed ecological niche models, among other names) provide a means to project changes in species' distributions in response to climate change Guisan et al., 2017). Although SDMs have projected climatically suitable habitat for a rare reptile found only in NZ, the tuatara Sphenodon punctatus, the sole surviving rhynchocepalian reptile (Jarvie et al., 2021), the vast majority of NZ reptile diversity has not been studied in this way (but see Bellis et al., 2020 who used SDMs to evaluate climatic suitability as a predictor of conservation translocation failure for two NZ lizards under current climate and Stenhouse et al., 2018 who combined an empirical model of temperature-dependent embryonic development with a mechanistic model of soil temperatures to examine range extent for a NZ lizard). Here, we build SDMs for NZ geckos and skinks explicitly considering spatial sampling bias in occurrence records, model over-fitting and multi-collinearity in climatic variables using current climate (representative of 1979-2013), before projecting the SDMs to a future time period (2070, average for 2061-2080). We then investigated how climatically suitable area varies with unlimited dispersal and no-dispersal scenarios. We next use Bayesian phylogenetic mixed models (BPMMs) to assess whether: (1) larger bodied species are more vulnerable than smaller bodied species, (2) nocturnal species are more vulnerable than diurnal species, (3) species at higher elevation are more vulnerable than species at lower elevation and (4) species at higher latitude are more vulnerable than species at lower latitude. Finally, we discuss the relationships of projected change in climatic suitability to these life-history traits and current distributions of NZ lizards, and the implications of our findings for conservation of NZ geckos and skinks under climate change.

| MATERIAL S AND ME THODS
Our method includes three distinct steps for assessing vulnerability of endemic NZ geckos and skinks to climate change: (1) infer time-calibrated phylogenetic trees, (2) use SDMs to predict climatic suitability under current and future climates and (3) run BPMMs to examine climate change vulnerability accounting for phylogeny while testing for body size, nocturnality, current elevation and current latitude. These three steps are summarised below; further technical details are provided in the Supplementary Information.

| Phylogeny
Phylogenetic trees were previously constructed using gene sequences for most NZ geckos (Nielsen et al., 2011) and skinks (Chapple et al. 2009). We supplemented existing sequences where necessary to include more recently described species, with two mitochondrial and two nuclear genes for geckos and five mitochondrial and one nuclear gene for skinks. We used the program BEAST (Drummond et al., 2012) to infer time-calibrated trees using fossil and biogeographical calibrations. We then extracted maximum clade credibility (MCC) trees containing all but three skink species and two gecko species that are currently recognised or hypothesised (See Supplementary Information for phylogenetic details, Table S1 for substitution models, and Figures S1 and S2 for the MCC phylogeny of geckos and skinks, respectively).

| Occurrence data
We used occurrence records in our SDMs for extant NZ gecko and skink species mainly from the NZ Department of Conservation Herpetofauna Database (Department of Conservation 2020). This database is the most comprehensive repository of occurrence records for endemic NZ lizards. After updating the taxonomy and cleaning the data to reflect the taxonomy as at 2019 of 43 gecko species recognised across seven genera (Dactylocnemis, Hoplodactylus, Mokopirirakau, Naultinus, Toropuku, Tukutuku and Woodworthia) and 61 skink species in one genus (Oligosoma) (Chapple et al., 2009;Hitchmough, Barr, et al., 2016;Hitchmough, Patterson, et al., 2016;Nielsen et al., 2011), we were left with a total of 5739 occurrence records, with 2387 gecko records and 3354 skink records (see Supplementary Information for steps, Table S2 for the number of raw and cleaned occurrence records for the modelled species and Figure S3 for the locations of occurrence records).

| Climatic variables
As predictors in our SDMs, we used six climate variables selected with biological justification to have ecological relevance for NZ lizards from the CHELSA dataset (Karger et al. 2017). The temperature and precipitation variables from the CHELSA dataset perform better for the spatial prediction of variables than other high-resolution datasets (Karger et al. 2017). The six predictor variables were as follows: annual mean temperature (bio1), maximum temperature of the warmest month (bio5), minimum temperature of the coldest month (bio6), mean temperature of driest quarter (bio9), mean temperature of wettest quarter (bio10) and precipitation of the driest quarter (bio17). The ecological relevance of the variables were as follows: (1) annual mean temperature and mean temperature of the warmest quarter likely to influence embryonic development, and thus limit where species can successfully reproduce; (2) mean temperature of the driest quarter and precipitation in the driest quarter likely influence susceptibility to desiccation in NZ lizards, which appear to have high rates of water loss, and thus limit distribution of the species (Hare & Cree, 2016)

| Modelling species distributions
We modelled the distribution of New Zealand lizards using the maximum entropy algorithm (Maxent, v 3.4.1; Phillips et al., 2017), a presence-background approach based on an inhomogeneous Poisson process that is designed to deal with presence-only data such as from occurrence record databases (Elith et al., 2011). We chose Maxent due to its superior performance over other SDM algorithms that use presence-only data like used in our study here (Elith et al., 2006), especially for species with small sample sizes (Pearson et al., 2006). To reduce the effects of sampling bias, we applied a weighted-target group approach (Phillips et al., 2009) using as background data any occurrence records for a NZ lizard ( Figure S3) that fell within 150 km of a presence record. This region includes environments that are probably accessible to NZ lizards given their dispersal limitations and configuration of barriers (Anderson & Raza, 2010;Barve et al., 2011). We balanced model fit with predictive ability by determining the optimal regularisation multiplier and feature class settings (Muscarella et al., 2014). Specifically, we built models with combinations of regularisation multiplier values from 0.5 to 4, in increments of 0.5, and the following feature-class combinations: linear; linear and quadratic; and linear, quadratic and product. We tested only these combinations of regularisation multipliers and feature classes to avoid overfitting and to produce interpretable response curves (Merow et al., 2013). Although ensemble models can reduce uncertainty caused by using a single SDM, recent research shows tuned models can perform comparably to ensembles (Hao et al., 2020). To reduce multicollinearity in the SDM for each species, we only retained climate variables with a variance inflation factor <10 (Dormann et al., 2013) within the 150 km buffer. Because the study involved model transfer across space and time, we used the 'block' method to partition localities (Roberts et al., 2017). From model outputs fitted with a clog-log transformation, we generated species distribution predictions for current and future climates. We then thresholded these probabilistic suitability maps to binary projections of suitable and unsuitable habitat for current and future climate scenarios using the maximum sum of sensitivity and specificity, a frequently recommended approach that tends to reflect the prevalence of modelled species well (Liu et al., 2016). For future climate scenarios, we generated consensus maps for each RCP using unweighted ensemble means from the respective GCMs.
To evaluate the performance of the optimal SDMs, we used the continuous Boyce index (CBI), and to assess model fit we used the test point omission rate based on the minimum training presence value (OR MTP ). The CBI represents the correlation between the model predictions and true probability presence (Boyce et al., 2002;Hirzel et al., 2006). Values of CBI range from −1 to 1, with values >0 indicating more accurate models and values <0 indicating performance worse than random. The threshold-dependent metric OR MTP indicates the proportion of testing localities with suitability values lower than that associated with the lowest-ranking training locality . The OR MTP values range from 0 for models that are not overfitted to 1 for models that are overfitted.

| Statistical analysis
For each species, we calculated the projected relative change in climatically suitable area under climate change. We also examined a nodispersal scenario by restricting future projections to areas identified as being currently climatically suitable. This no-dispersal scenario is relevant because many NZ lizards do not disperse large distances or over highly fragmented habitat . We then used the proportion of pixels present from the binary SDMs to build potential scenarios (e.g., future projected range or the number of cells projected to be occupied/current predicted range or the number of cells predicted to be occupied), with dispersal and no-dispersal scenarios. As a metric of climate vulnerability, we calculated a log range ratio (hereafter lnRR), the natural log-transformed ratio of the projected range area in 2070 to the present range area (after adding 1 to both the numerator and the denominator, both of which are in units of km 2 ).
We investigated projected changes in climatically suitable area with BPMMs to account for non-independence in species' characteristics due to common ancestry. We carried out the BPMMs with MCMC estimation using the R package MCMCglmm (Hadfield, 2010) for species with matching phylogenetic data and SDMs. We used separate BPMMs with a Gaussian error structure with identity link for both geckos and skinks to evaluate whether lnRR is related to phylogeny, distribution or life-history traits. We incorporated four fixed effects in our BPMMs: body size (continuous variable), nocturnality (two-level categorical variable: diurnal and nocturnal), current mean latitude (continuous variable) and current mean elevation (continuous variable). Body size and activity phase information were compiled from the Landcare Research New Zealand Lizards Database (Bell, 2010), with body size not separated by sex as there is no substantial sexual dimorphism in NZ lizards (Chapple, 2016). The body size of lizard species ranged in snout-vent length (SVL) from 51 mm (O. levidensum) to 161 mm (H. duvaucelii). Although activity phase of lizards has often been classified dichotomously (diurnal v nocturnal, sometimes crepuscular), these terms do not always do justice to the activity patterns of NZ lizards (reviewed in Hare and Cree (2016)).
Nevertheless, as a simple descriptor activity phase, that is, the main activity/foraging times, can still be a useful predictor when assessing the vulnerability of reptiles to climate change (e.g., Meiri et al., 2013).
Current mean latitude was estimated from the occurrence records and current mean elevation was extracted at the locations of the occurrence records from the ca. 90 m EarthEnv-DEM90 digital elevation model (Robinson et al., 2014). Mean latitude and mean elevation were included to establish whether species in certain regions were more susceptible to climate change as has been shown in previous studies (e.g., Parmesan, 2006). To improve interpretation of regression coefficients, we standardised the continuous input variables by subtracting the mean and dividing by two standard deviations (Gelman, 2008). We report the estimates of all regression coefficients (β) as posterior means with their 95% credible intervals (CIs), and consider effects statistically significant when the 95% CIs did not overlap zero (see Supplementary Material for further details on BPMMs.) We included in our BPMMs output modelled distributions for NZ lizards where we ran SDMs with greater than 15 spatially thinned occurrence records as Maxent can accurately predict distributions for species with small sample sizes (Pearson et al., 2006), although we also ran a sensitivity analysis where we excluded species with fewer than 30 spatially thinned occurrence records as below this threshold the performance of SDMs can be reduced for model projections (Wisz et al., 2008). As the results from our sensitivity analyses were qualitatively similar, we report findings from SDMs trained with >15 filtered occurrence records in the main text (see Supplementary Material for results from sensitivity analyses).

| Phylogeny reconstruction
Time-calibrated phylogenies estimated with BEAST were highly consistent with previous reconstructions for both NZ geckos and skinks ( Figure S2 for geckos and Figure S3 for skinks).
For skinks, the area of climatic suitability decreased for 29 species with RCP4.5 and for 25 species with RCP8.5 (Table S3; Figure 2a,b).

| Phylogenetic mixed-effect models
All BPMMs performed better when phylogeny was included (∆DIC > 2). For geckos, the projected relative changes in climatic suitability area (i.e., lnRR) showed non-significant relationships for life-history traits (body size and activity phase) and latitude.

F I G U R E 1
Maps of species richness for New Zealand diplodactylid geckos and eugongylinae skinks generated using maximum entropy (Maxent) models under current (1979-2013) and future climate scenarios (2070, average 2061-2080) for Representative Concentration Pathways (RCPs) 4.5 and 8.5. We also show a no-dispersal scenario for geckos and skinks under RCP4.5 and RCP8.5 However, relationships with elevation were positively significant (Figure 3a,b). For the no-dispersal scenario, life-history traits and current distribution showed similar relationships (Figure 3a,b). The relationships were the similar for both emission scenarios, although there was more variation for the high emission scenario of RCP8.5 than the low emission scenario of RCP4.5.
For skinks, the projected relative change in climatic suitability area (i.e., lnRR) showed non-significant relationships with activity phase, elevation and latitude but a negative significant relationship for body size (Figure 3a,b). The no-dispersal showed similar relationships (Figure 3a,b). Relationships were similar for both RCP4.5 and 8.5, with more variation for the high emission scenario.

| DISCUSS ION
Human-mediated climate change threatens species with declines in climatically suitable area and possible extinction (Pecl et al., 2017;Urban, 2015). Our SDM projections showed variable vulnerability to climate change for NZ lizards, with climatically suitable area F I G U R E 2 Percentage of area change in climatic suitability for New Zealand diplodactylid geckos and eugongyline skinks as generated using maximum entropy (Maxent) models under current  and future climate scenarios (2070, average 2061-2080) for Representative Concentration Pathways (RCPs) 4.5 and 8.5 (panels a and b, respectively). We also show a nodispersal scenario for geckos and skinks under RCP4.5 and RCP8.5 (panels c and d, respectively) F I G U R E 3 Regression coefficients from Bayesian phylogenetic mixed models (BPMMs) for New Zealand diplodactylid geckos (panels a and b) and eugongyline skinks (panels c and d), showing both a relative change in climatically suitability habitat and a no-dispersal scenario under current  and two future climates (2070, average 2061-2080) via the consensus of four general circulation models for the Representative Concentration Pathways (RCPs), RCP4.5 and RCP8.5. As the response variable in the BPMMs, we used log-transformed ratio of the projected range area relative to the current range area (lnRR for log range ratio) while the predictors variables were body size, activity phase (diurnal or nocturnal), mean elevation and mean latitude. Mean estimates from regression coefficients are shown by circles for relative change in climatic suitability (black) and no dispersal (grey) with 95% credible intervals shown by lines. See Figure S4 for regression coefficients from BPMMs trained with the results of the SDMs for NZ lizards with ≥30 spatially thinned occurrence records decreasing for most, but not all, geckos and skinks. Our BPMM results showed phylogenetic signal in the projected change in climatically suitable area for NZ lizards, suggesting more closely related species will be more similarly impacted by climate change.
Additionally, geckos at higher elevation were predicted to be more impacted, as were smaller bodied skinks. Our no-dispersal scenario showed qualitatively similar results for both the phylogenetic signal and the predictor variables, although the amount of climatic suitability area substantially decreased relative to a scenario allowing for dispersal. Below, we discuss the projected change in climatic suitability identified by the SDMs for NZ lizards, the relationships of projected change in climatic suitability to life-history traits and current distributions, and the implications of our findings for conservation of NZ lizards under climate change.

| Species distribution models
Our findings show variable vulnerability to projected change in climatic suitability among NZ lizards under climate change. Although the amount of climatically suitable area increased for a minority of NZ lizards under climate change, it decreased for the majority and disappeared entirely under the business-as-usual scenario for three species, the geckos W. 'Otago Southland large' and T. rakiurae and the skink O. homalonotum (Figures 1 and 2b). Our results showed that climatic suitability was projected to generally shift poleward for NZ lizards (Figure 1), a finding similar to a meta-analysis showing that, on average, the ranges of terrestrial organisms have moved poleward in response to ongoing climatic shifts (Chen et al., 2011). However, there was also some variation in the direction of change in climatic suitability from the projections of the SDMs for NZ lizards, as has been documented empirically for other species due to climate-related range shifts (Lenoir & Svenning, 2015). For the no-dispersal scenario, projections of climatic suitability for most species did not overlap with their current range (Figures 1, 2c,d), indicating that without dispersal through areas in which they have not been detected to date, many species are more vulnerable to climate change. A recent global-scale study similarly showed that dispersal ability for a large number of terrestrial species across taxonomic groups (invertebrates, vertebrate and plants) was critical to prevent widespread losses of climatic suitability area in the future (Warren et al., 2018). A recent regional study of mammals in Ecuador also showed the importance of dispersal in predicting changes in species diversity, with inability to disperse substantially impacting climatically suitable area (Iturralde-Pólit et al., 2017).
We found that climatically suitable area decreased considerably if we followed the current, business-as-usual trajectory for RCP-8.5, compared to the stabilisation strategy for RCP-4.5 (Figure 2a,b).
This finding supports earlier studies suggesting that future global extinction risk and population declines of a wide range of taxa are predicted to increase as global temperatures rise (Urban, 2015;Warren et al., 2018). For NZ lizards, when climatic space for a species is projected to decline in areas where they are currently found, conservation translocations-the intentional movement and release of organisms to restore populations-could be used to establish populations in regions that will remain climatically suitable into the future (IUCN/SSC, 2013). Conservation translocations could include reintroductions, the re-establishment of the focal species within its indigenous range, and managed relocations, the movement of the focal species outside its indigenous range to avoid extinctions of populations . Due to the vulnerability of many NZ lizard species to predation from introduced mammalian predators, priority areas for conservation translocations should be to destination sites where the full suite of invasive mammalian predators are controlled to very low levels or have been eradicated (Nelson et al., 2014;Towns et al., 2016). Opportunities for conservation translocations of NZ lizards have focused on reintroducing species to offshore islands as their introduced mammalian predators have been eradicated (Nelson et al., 2014) and, more recently, to predatorfree reserves on the mainland (i.e., North and South Islands) where these predators have been eradicated within mammal-proof fences (Jarvie et al., 2014(Jarvie et al., , 2015Nelson et al., 2014;Reardon et al., 2012). Where climatically suitable area is projected to disappear entirely, or to decline substantially (Figure 2), individuals from species could be bought into captivity to safeguard against the potential threat of future climate change while further research is undertaken into climate tolerance (e.g., Chukwuka et al., 2020). Populations in captivity could function as safe havens for vulnerable species threatened by climate change, as conservation biologists and natural resource managers seek to better understand the impact of ongoing global environmental change on wildlife (Conway, 2011). Risks assessment for potential conservation translocations, particularly for managed relocations, should be undertaken to navigate the complex ethical challenges arising from novel interactions (Burbridge et al., 2011) and risks of collateral damage (Ricciardi & Simberloff, 2009).
Despite SDMs for some NZ lizards being trained with a relatively small number of filtered occurrence records, that is, for 11 of the 51 modelled species with >15 records but <30 records, the number used is commonly understood to be sufficient to allow robust analysis (Pearson et al., 2006;Warren et al., 2018), as well as an algorithm in Maxent whose performance is less sensitive to sample size (Wisz et al., 2008). Furthermore, our sensitivity analyses where we compared BPMMs from the output of SDMs trained with ≥30 filtered occurrence records with the SDMs with >15 filtered occurrence records produced qualitatively similar results. We also accounted for the cryptic nature of many NZ lizards (Chapple, 2016) through the use of a weighted-target group approach (Phillips et al., 2009). This widely used approach deals with sampling bias for presence-only SDMs better than traditional approaches with random sampling (Syfert et al., 2013).
Our choice with the SDMs to use regularisation multipliers and feature classes to avoid overfitting (Merow et al., 2013) as well as the 'block' method to partition localities should improve performance from the projections of climatic suitability across space and time (Roberts et al., 2017). Although we acknowledge that some species have already experienced substantial anthropogenic range contractions (Nelson et al., 2014), only current records are available for most NZ lizards.
Anthropogenic range contractions have been accounted for in SDMs where a detailed vertebrate fossil record exists (e.g., Faurby & Araújo, 2018;Jarvie et al., 2021), or with estimated range maps (e.g., Jarvie & Svenning, 2018;Monsarrat et al., 2019). In the NZ fossil record, smaller squamates have not usually been identified to the species level (Worthy, 2016), making this approach difficult for the majority of NZ lizards. The use of ancient DNA with fossil records from NZ lizards could also provide a much better understanding of past geographical ranges for smaller squamates (e.g., Seersholm et al., 2018). A promising alternative method to the use of the above approaches in correlative SDM to overcome anthropogenic range contractions would be the use of mechanistic models to forecast potential distributions. Mechanistic models often integrate physiological tolerances and/or demographic parameters of species to climate (Kearney & Porter, 2009); however, they require intensive parameterisation and therefore they are to date only available for a few well-studied species (Briscoe et al., 2019).

| Phylogenetic mixed-effect models
Our BPMM analyses showed for geckos that the projected relative changes in climatic suitability area (i.e., lnRR) were significantly positively related to elevation but not significantly related for the other predictor variables (Figure 3a Previous studies have shown that increased body size is associated with increased extinction risk (Dirzo et al., 2014), with a study of NZ lizards also suggesting heightened extinction risk associated with increased body size . This is due to larger bodied lizard species being more vulnerable to predation by introduced mammalian predators than smaller bodied species as the refugia they use when inactive are larger and therefore more accessible to mammalian predators (Nelson et al., 2014;Whitaker, 1978). A possible reason why body size in our analyses was a non-significant predictor for geckos, but was negatively significant for skinks was due to body size differing less for the geckos than skinks (Table S2). This could be because many, but not all, larger NZ skinks are nocturnal, whereas smaller species are predominantly diurnal.
The larger bodied species, both geckos and skinks, have suffered more historic range contractions than smaller ones (Nelson et al., 2014).
In a recent review, high-elevation locations were noted as being among the most sensitive ecosystems to climate change, with species found in such locations being affected at a faster rate than other terrestrial habitats (Pecl et al., 2017). Our findings showing geckos being impacted by climate change at higher elevations, but not skinks, could be because we were able to run SDMs for more alpine geckos that live at higher elevations relative to skinks (due to low number of occurrence records for many alpine skink species).
Further research is needed to explore the risk of climate change to NZ's alpine lizards (but see, e.g., Chukwuka et al., 2020).
We found phylogenetic signal in the BPMMs for projected change in climatic suitability for NZ lizards, implying that closely related species are more likely to show similar trends. A previous study of plant, bird and mammal assemblages across Europe using a consensus of ensemble forecasts and high-resolution phylogenetic trees showed the consequence of climate change on phylogenetic diversities, with species vulnerability to climate change weakly clustering across phylogenies (Thuiller et al., 2011). Although projections in this previous study did not predict a large drop in phylogenetic diversity, they did suggest a future restructuring of the spatial distribution of the tree of life (Thuiller et al., 2011). Our findings similarly suggest phylogenetic diversity for NZ lizards will likely be restructured under climate change, but the extent largely depends on species' ability to disperse ( Figure 1). Future studies should investigate in more detail for NZ lizards how both species richness and phylogenetic diversity could shift under climate change (e.g. Di Virgilio et al., 2017), with results being used to inform conservation management of threatened species.

| Conservation management of NZ lizards under climate change
Predictive distribution models are increasingly used to support conservation decision-making (Guisan et al., 2013), including identifying critical habitats and mapping suitable sites for conservation translocations. Our results suggest that the majority of NZ lizards will be negatively impacted by the direct effects of climate change, with species that were projected to experience declines showing reductions in climatically suitable area of 42%-46% for geckos and of 33%-52% for skinks ( Figure 2). The declines in climatically suitable area only increased further for species when we modelled a no-dispersal scenario. Note that although we modelled climatic suitability for many NZ lizards (51 of the 104 species, with 28 of the 61 skinks and 23 of the 43 geckos), it is likely that many of the species we did not run SDMs for (due to the limited number of occurrence records) will be especially vulnerable to climate change, as they are occurring at a single location or are range restricted (<1000 km 2 ; Townsend et al., 2008;Hitchmough et al., 2021). Such species should be closely monitored to assess the potential adverse effects from climate change, including continued investigations into the possible heightened vulnerability for species with large body size and primarily nocturnal activity phase. The use of conservation translocations should also be considered if there is a risk their naturally restricted range will become unsuitable. For the modelled NZ lizards with projected declines in climatic suitability, we recommend translocations following further assessments to identify climatically suitable areas, more detailed research into climate tolerance, andparticularly if the two former actions were deemed unfeasible-that individuals should be brought into captivity to establish ex situ populations. Such conservation interventions should help to safeguard NZ's endemic lizard fauna into the future.

| CON CLUS ION
Our findings that NZ lizards vary in vulnerability to climate change can be used to inform conservation management for species. The results can be used to identify climatically suitable areas for species under current and future climates by prioritising in situ conservation, identifying regions suitable for conservation translocations and identifying ex situ conservation opportunities. More broadly, our research highlights the need to assess climatic suitability for vulnerable ectotherm faunas under future climate scenarios. Lizards in the temperate regions of the Southern Hemisphere can be highly speciose but understudied relative to tropical and Northern Hemisphere lizards. The high endemicity and relatively small land mass available to southern temperate lizards make it vital that we understand how they will be affected by ongoing global environmental change.

ACK N OWLED G EM ENTS
For occurrence records and tissue samples, we thank Te Papa Atawhai | Department of Conservation (DOC), Tony Jewell, Marieke Lettink, Dylan van Winkel and numerous others. For this particular research, collecting permits were not needed. We thank Alison Cree for early discussions on the use of species distribution models, and DOC, regional councils and iwi (Māori tribes) for support. We also thank the editorial team, including Christine Meynard, Ge-Xia Qiao and Michael Dawson, and two anonymous reviewers for providing constructive comments on earlier versions of the manuscript.
Funding was provided by a DOC contract to SJ. DGC was supported by a grant from the Australian Research Council (FT200100108)

DATA AVA I L A B I L I T Y S TAT E M E N T
Scripts to run examples of species distribution models and Bayesian phylogenetic mixed models in the programming language R are provided in the Supporting Information. We also provide MCC tree files for the geckos and skinks. Although GPS coordinates for current populations are not included due to the potential threat of poaching, the climate variables for each species are provided in the Supporting Information. The climate data were from the CHELSA database (https://chels a-clima te.org/), which can be freely downloaded for current and future scenarios.