Assessment of the vulnerability of Spiny Lizards in Central Mexico under anthropogenic pressure: integrating habitat suitability, landscape connectivity, fragmentation and protected area coverage

351 Assessment of the vulnerability of Spiny Lizards in Central Mexico under anthropogenic pressure: integrating habitat suitability, landscape connectivity, fragmentation and protected area coverage Jossye Vargas-Jaimes1, Javier Manjarrez2, Fabiola Ramírez-Corona3, Guillermo A. Woolrich-Piña4, María G. González-Pedroza2, María Elena Estrada-Zúñiga2, Carlos Alejandro Rangel-Patiño5, Crystian Sadiel Venegas-Barrera6, Aurelio Nieto-Trujillo2, Alicia Monserrat Vazquez-Marquez2, Juan Carlos Guido-Patiño1, Armando Sunny1 1 Centro de Investigación en Ciencias Biológicas Aplicadas, Facultad de Ciencias, Universidad Autónoma del Estado de México, Instituto literario 100, Colonia Centro, 50000, Toluca, Estado de México, Mexico 2 Facultad de Ciencias, Universidad Autónoma del Estado de México, Instituto literario 100, Colonia Centro, 50000, Toluca, Estado de México, Mexico 3 Taller de Biogeografía y Sistemática, Facultad de Ciencias, Universidad Nacional Autónoma de México, Av. Ciudad Universitaria 3000, C.P. 04510, Coyoacán, CDMX, Mexico 4 Laboratorio de Zoología, División de Biología, Subdirección de Investigación y Posgrado, Instituto Tecnológico Superior de Zacapoaxtla, Carretera AcuacoZacapoaxtla km. 8, Col. Totoltepec, C. P. 73680, Zacapoaxtla, Puebla, Mexico 5 Tecnológico de Estudios Superiores de Huixquilucan, Paraje El Río S/N, La Magdalena Chichicaspa, Huixquilucan, 52773, Mexico 6 Instituto Tecnológico de Ciudad Victoria, Tecnológico Nacional de México, Ciudad Victoria P.C. 87010, Mexico Corresponding author: Armando Sunny (sunny[email protected]) Copyright: © Jossye Vargas-Jaimes et al. This is an open access article distributed under terms of the Creative Commons Attribution License (Attribution 4.0 International – CC BY 4.0). Research Article Abstract The genus Sceloporus represents a diverse and ecologically significant group of lizards within the Mexican herpetofauna, particularly in the Trans-Mexican Volcanic Belt (TMVB), a region of exceptional endemism increasingly threatened by anthropogenic pressures. In this study, we conducted a comprehensive spatial vulnerability assessment of 17 species of Sceloporus that are present in the TMVB, integrating ecological niche modelling (ENM), habitat fragmentation metrics, land-use exposure and protected area overlap. Binary distribution maps were generated from the ensemble models optimised with biomod2. For each species, we calculated the habitat area, number of habitat patches and mean patch size, synthesising these values into a composite vulnerability index. Land-use pressures were quantified through spatial overlaying of urban, agricultural and pasture-land layers, followed by principal component analysis (PCA) and k-means clustering to generate a multivariate anthropogenic pressure zonation. Species distributions were also compared against those of state and federal protected areas. The 17-species analysis revealed noteworthy variation in spatial dimensions and fragmentation of habitat, exhibiting high spatial vulnerability and high exposure to anthropogenic pressure. Furthermore, our results revealed that the predicted habitats of all 17 species overlap with highly transformed urban landscapes. A bivariate risk assessment identified S. megalepidurus as the only species with both high pressure and high vulnerability, while the species S. anahuacus, S. albiventris and S. minor, currently listed as “least concern” or “not evaluated” by the IUCN, showed alarming risk signals when spatially explicit metrics were incorporated. Our integrative framework underscores the need to revisit conservation assessments by incorporating habitat quality, fragmentation and exposure to land-use pressures, not just geographic range. Academic editor: Clarissa Rosa Received: 20 May 2025 Accepted: 6 September 2025 Published: 3 October 2025 ZooBank: https://zoobank.org/ C3582276-E64E-4616-9B2753D1194E0F10 Nature Conservation 59: 351–386 (2025) DOI: 10.3897/natureconservation.59.159660 352 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards Key words: Anthropogenic pressure, conservation prioritisation, ecological niche modelling, Sceloporus species, Trans-Mexican Volcanic Belt, vulnerability index Introduction Mexico ranks second worldwide in reptile diversity, with 1,031 documented reptile species, nearly half of which are endemic (Cabrera-Hernández et al. 2024; Uetz 2025). Despite this remarkable richness, many species face significant conservation challenges, including those classified as Data Deficient or Not Evaluated on the IUCN Red List – categories that inherently increase their vulnerability (Suazo-Ortuño et al. 2023). These risks are further compounded by factors such as restricted distribution ranges, dependence on specialised microhabitats and intensifying anthropogenic pressures (Rosas-Espinoza et al. 2024). Such threats are particularly pronounced in the Trans-Mexican Volcanic Belt (TMVB), a region where biodiversity and human activity intersect. The TMVB spans Mexico from west to east and is widely recognised as one of the country’s most ecologically complex and biogeographically significant regions (Smith 1941; Bryson et al. 2011; Gámez et al. 2012). Its diverse topography, pronounced climatic gradients and dynamic geological history have fostered exceptional levels of species richness and endemism, especially amongst amphibians and reptiles (Flores-Villela and Canseco-Márquez 2007; Sunny et al. 2017). Indeed, the TMVB ranks second in overall herpetological richness in Mexico and holds the highest number of endemic reptile and amphibian species (Flores-Villela and Canseco-Márquez 2007). Approximately a quarter of the nation’s herpetofauna, more than 80% of which are endemic, are found within this biogeographic region (Lemos-Espinal and Smith 2024), underscoring its crucial role in the conservation of Mexico’s vertebrate biodiversity. However, the ecological integrity of the TMVB is increasingly threatened by expanding urbanisation, agricultural intensification, deforestation and habitat fragmentation (Ochoa-Ochoa and Flores-Villela 2006; CONAPO 2010; Sunny et al. 2017; INEGI 2020). These pressures are exacerbated by the fact that the TMVB encompasses some of the most densely populated and industrially active regions in Mexico (González-Fernández et al. 2018). Although the establishment of protected natural areas (PNAs) has aimed to mitigate biodiversity loss, many of these areas are insufficiently managed or fail to represent key ecosystems and species (Ruiz-Reyes et al. 2024; Venegas-Barrera et al. 2024). As a result, large portions of biologically rich landscapes remain unprotected and are undergoing rapid environmental change (Sinervo et al. 2010; González-Fernández et al. 2022; Sunny et al. 2024a, 2024b). Sinervo et al. (2010) demonstrated that climate change has already driven substantial local extinctions in lizard populations, with 12% of Mexican populations disappearing since 1975. Their work validated physiological extinction risk models with empirical field data, projecting that, by 2080, local extinctions could reach 39% globally, with up to 20% of species lost. These alarming trends, confirmed across multiple continents, indicate that lizards have already crossed a threshold of climate-induced extinctions. When compounded with habitat degradation, fragmentation and loss of ecological connectivity, these climate-driven pressures threaten to accelerate biodiversity loss, particularly in ecologically sensitive and poorly-protected Citation: Vargas-Jaimes J, Manjarrez J, Ramírez-Corona F, Woolrich-Piña GA, GonzálezPedroza MG, Estrada-Zúñiga ME, Rangel-Patiño CA, VenegasBarrera CS, Nieto-Trujillo A, Vazquez-Marquez AM, GuidoPatiño JC, Sunny A (2025) Assessment of the vulnerability of Spiny Lizards in Central Mexico under anthropogenic pressure: integrating habitat suitability, landscape connectivity, fragmentation and protected area coverage. Nature Conservation 59: 351–386. https://doi.org/10.3897/ natureconservation.59.159660 353 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards regions, such as the Trans-Mexican Volcanic Belt (Mastretta‐Yanes et al. 2015; González-Fernández et al. 2018; Sunny et al. 2022; Rubio-Blanco et al. 2024). Understanding how species respond to these pressures is essential for effective conservation planning (Tischendorf and Fahrig 2000; Bennett 2003; Guisan and Thuiller 2005; Guisan et al. 2013; Maciel-Mata et al. 2015; Ersoy et al. 2018). This is particularly important for taxa, such as the genus Sceloporus, which is emblematic of Mexican herpetofauna. With approximately 110 recognised species – 106 of which are endemic to Mexico (Wiens et al. 2010; Uetz 2025) – Sceloporus represents a highly diverse lineage of spiny lizards that occupy a wide range of habitats, including montane forests, arid scrublands and grasslands. Despite the ecological importance and diversity of this group, there are significant knowledge gaps regarding the environmental factors that influence their distribution, as well as their vulnerability to anthropogenic pressures (Smith 1941; Sites et al. 1992; Wiens and Reeder 1997; Smith and Ballinger 2001; Flores-Villela and Canseco-Márquez 2004; Leaché 2010; Leaché et al. 2016; Martínez-Méndez et al. 2019). Moreover, although efforts have been made to characterise the thermal physiology, systematics and biogeographic history of Sceloporus species, fewer studies have integrated ecological niche modelling with multivariate pressure mapping, landscape fragmentation analysis and protected area overlap assessments (Leaché et al. 2010, 2016; López-Reyes et al. 2024; García-Rosales et al. 2024). These integrative approaches are necessary to evaluate whether current conservation instruments adequately protect the most vulnerable species and to identify areas where conservation gaps persist. The TMVB, with its high levels of endemism and strong anthropogenic pressure, is an ideal natural laboratory for implementing these tools and advancing conservation science. Despite their broad distribution, most Sceloporus species included in this study are currently classified as “Least Concern” (LC) by the IUCN Red List or have not been evaluated. Amongst them, S. megalepidurus is listed as Vulnerable (VU), while S. anahuacus, S. bicanthalis, S. bulleri, S. dugesii, S. grammicus, S. horridus, S. mucronatus, S. palaciosi, S. scalaris and S. torquatus are considered LC. In contrast, S. subniger, S. albiventris, S. melanogaster, S. minor, S. utiformis and S. variabilis have not yet been assessed (IUCN 2025). According to the Mexican official normative (NOM-059-SEMARNAT 2025), only S. grammicus and S. megalepidurus are listed as species under special protection (SEMARNAT 2025). This disparity underscores a critical limitation of global conservation assessments, which often rely on broad-scale distribution data and may underestimate species’ actual vulnerability to localised anthropogenic pressures. Addressing these gaps requires regional, fine-scale evaluations that integrate multiple indicators, such as ecological niche modelling, habitat fragmentation, land-use exposure and representation within protected areas. In this context, our study assessed the vulnerability of 17 Sceloporus species occurring in the Trans-Mexican Volcanic Belt (TMVB), providing updated insights into their conservation status. Our main objectives were to: (1) model the potential distribution of each species using ecological niche modelling techniques, including both Maxent via the ENMeval package and an ensemble approach integrating eight algorithms implemented in biomod2; (2) quantify the degree of habitat fragmentation and patch isolation using landscape metrics; (3) assess the overlap between species distributions 354 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards and anthropogenic land-use pressures, including urban areas, agriculture and cultivated pastures; and (4) evaluate the degree to which each species is represented within federal and state protected areas. We hypothesised that species with smaller and more fragmented distributions would be more exposed to anthropogenic land use and have lower protection coverage, thereby facing a compounded risk of extinction. Materials and methods The Trans-Mexican Volcanic Belt (TMVB) is a Neogene volcanic arc covering approximately 160,000 km2 across central Mexico that extends from 18°30'N to 21°30'N (Ferrari et al. 2018). Formed by the subduction of the Cocos and Rivera plates beneath the North American Plate, the TMVB is characterised by a chain of stratovolcanoes, volcanic fields and extensive highlands (Gómez-Tuena et al. 2005; Ferrari et al. 2012). Its rugged topography reaches elevations exceeding 4,000 m above sea level (m a.s.l.). The region’s complex geological history – marked by multiple volcanic episodes, tectonic activity and climatic fluctuations – has shaped a mosaic of habitats ranging from temperate forests and alpine grasslands to xerophytic scrublands (Mastretta-Yanes et al. 2015). These environmental gradients have fostered exceptional levels of species richness and endemism, making the TMVB a key biogeographic corridor and a biodiversity hotspot of global significance (Flores-Villela and Canseco-Márquez 2007; Morrone 2014). Despite its importance, several aspects of its geological evolution, ecological connectivity and response to anthropogenic pressures remain poorly understood. To investigate patterns of habitat suitability and conservation within this region, we compiled occurrence data for 17 lizard species of the genus Sceloporus, each with distributions overlapping portions of the TMVB (Arriaga et al. 1997). Species selection followed the criteria proposed by Campillo-García et al. (2021), particularly regarding taxonomic updates within the S. torquatus complex. Occurrence records were sourced from the iNaturalist online repository (https://www.inaturalist.org/) and consolidated from multiple biodiversity databases, including REMIB, UNIBIO, GBIF, HERPNET, VertNet and IREKANI, ensuring a comprehensive and updated dataset. To ensure data currency, only records from 2000 to 2023 were considered. To ensure model accuracy and mitigate potential sampling biases that might induce model overfitting (Boria et al. 2014), the obtained data were filtered using the NicheToolBox online tool (Osorio-Olvera et al. 2020). Duplicate observations within the same pixel, approximately 1 km2 in resolution, were eliminated, employing a δ threshold value of 0.0083333. Furthermore, records lacking precise coordinates, sampling year information or data falling beyond Mexican territorial confines were systematically excluded. After filtering, the dataset included, 61 for S. albiventris, 34 for S. anahuacus, 44 for S. bicanthalis, 17 for S. bulleri, 63 for S. dugesii, 621 for S. grammicus, 238 for S. horridus, 40 for S. megalepidurus, 132 for S. melanogaster, 235 for S. minor, 140 for S. mucronatus, 22 for S. palaciosi, 97 for S. scalaris, 122 for S. subniger, 255 for S. torquatus, 111 for S. utiformis and 440 for S. variabilis (See Suppl. material 1: table S1). The dataset is publicly available at: Zenodo. https://doi.org/10.5281/zenodo.16879333. 355 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards Environmental variables Three categories of environmental variables were employed to construct the ecological niche models: land use and vegetation, topographic and climatic. Land-use and vegetation data were sourced from the National Institute of Statistics and Geography for the year 2016 (INEGI 2016), while topographic data were obtained from a digital elevation model at a scale of 1:50000 (INEGI 2016). Climate data encompassing 19 climatic and 12 solar radiation variables were sourced from WorldClim (Hijmans et al. 2005), with a spatial resolution of 30 arc-seconds (approximately 1 km2). The land-use map underwent reclassification to delineate different classes, which were subsequently transformed into continuous variables using the filter module of IDRISI SELVA 17.0 (Clark Labs 2012). All layers were processed in raster format with a resolution of 1 km2 utilising the raster (Hijmans and Van Etten 2016) and rgdal (Bivand et al. 2015) packages in R software (version 4.4.2; R. Development Core Team 2024). Furthermore, Pearson’s correlation analysis was conducted to identify and eliminate highly correlated variables (r2 > 0.7; Dormann et al. (2013)), facilitated by the ENMtools and usdm packages in EcoNicheS (Sunny et al. 2025) for R software. The resulting layers were utilised in the preparation of ecological niche models (ENMs; See Suppl. material 1: table S2). Accessibility areas We generated the accessibility area (M) for each species using ArcMap software and the R programming language (ESRI 2020; R Studio Team 2021). We used delimitation to follow the biogeographic region containing the records of each species in Mexico (Arriaga et al. 1997; Sunny et al. 2017, 2024a, 2024b; González-Fernández et al. 2018; Rubio-Blanco et al. 2024) (Fig. 1, See Suppl. material 1: fig. S1). Ecological niche modelling To ensure a robust and comparative framework for ecological niche modelling, we employed two complementary approaches: Maxent modelling through the ENMeval package and ensemble modelling using the biomod2 package. The use of Maxent via ENMeval allowed us to explore the influence of different parameter settings (feature classes and regularisation multipliers) on model performance and spatial predictions, providing a detailed understanding of how algorithm configuration can affect outcomes. Simultaneously, we implemented ensemble models in biomod2, integrating eight commonly used algorithms (including Maxent) to generate consensus predictions that are generally more stable and less sensitive to the assumptions of individual methods. By comparing the results obtained from these two modelling strategies, we aimed to assess the consistency of spatial patterns across different methodological frameworks. This dual approach not only strengthens the reliability of our findings, but also highlights potential areas of uncertainty, which is critical when providing information for conservation strategies, based on species distribution models. 356 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards Figure 1. Accessibility areas (M) for each species (dark grey polygons) and occurrence records used in ecological niche models for 17 species of Sceloporus. Coloured dots represent occurrence data, with colours corresponding to the spatial partitions used for model calibration and evaluation. For most species, the block partitioning method was applied; however, for S. bulleri and S. palaciosi, a jackknife partitioning method was used due to the low number of records. In the upper right corner, a map of Mexico (in grey) shows the study area M for each species, with the Trans-Mexican Volcanic Belt (TMVB) highlighted in green. 357 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards ENMeval: Maximum entropy algorithm The maximum entropy algorithm was implemented using Maxent 3.3.3 software (Phillips et al. 2006) in conjunction with species-specific optimisations facilitated by the ENMeval package (Muscarella et al. 2014) in R software (R Studio Team 2021). The construction of species-specific ecological niche models was based on occurrence records and preselected environmental variables. To ascertain the optimal model settings, a preliminary test was conducted using the ENMeval package, encompassing all feasible combinations of parameters. This pilot test selected regularisation values of 1 and 2, along with characteristic classes, including “linear” (L), “quadratic” (Q), “product” (P), “threshold” (T) and “hinge” (H), encompassing combinations, such as L, H, LQ, LQP, LQH, LQHP and LQHPT. Additionally, 100,000 random points were generated as background points during the modelling process by utilising partitioning methods, namely, “Block”, “Checkerboard 1” and “Checkerboard 2”. For species with limited records, such as S. bulleri and S. palaciosi, the “jackknife” data partitioning method was employed. Optimal parameter combinations for each species were determined, based on the corrected Akaike Information Criterion (AICc) (Warren et al. 2010; Warren and Seifert 2011). Model accuracy and reliability were assessed using various approaches, including the area under the curve (AUC) and transformation of the AUC into a partial receiver operating characteristic (ROC) curve, facilitated by the ntbox package in R software. biomod2 An ensemble model was also constructed with the biomod2 package (Thuiller et al. 2025) implemented in the shinydashboard package EcoNicheS for R. Eight of the ten algorithms available in the package were selected. These models included the generalised linear model (GLM), generalised boosting model (GBM), generalised additive model (CTA), artificial neural network (ANN), surface range envelope or BIOCLIM (SRE), flexible discriminant analysis (FDA), random forest (RF), Maxent and MaxNet. As presence data were available, a set of 100000 pseudoabsence data points was randomly generated in EcoNicheS for R. The prevalence was set to 0.5 to ensure the presence or absence of the same importance in the calibration process (Barbet‐Massin et al. 2012; Linero et al. 2020). The models were calculated ten times, each run with a different selection of training and testing datasets. The proportion of data for calibration was set to 70% (Linero et al. 2020). To evaluate the overall model accuracy of both the Maxent and ensemble models, the area under the curve (AUC) from the receiver operating characteristic (ROC) plot was used first (Metz 1978; Phillips et al. 2006), followed by the transformation of the area under the curve (AUC) into a partial receiver operating characteristic (PROC) curve (partial ROC; Peterson et al. (2008)). For this purpose, the package EcoNicheS (Sunny et al. 2025) in R was used. For the ensemble model, the true skill statistic (TSS) was calculated with the package biomod2 (Allouche et al. 2006), based on the sum of sensitivity and specificity minus one (the sensitivity is the proportion of correctly predicted presences and the specificity is the proportion of correctly predicted absences). In the biomod2 approach, the models were assembled by a total consensus rule and algorithms were assembled by the weighted sum of evaluations greater than TSS > 0.8 (Linero et al. 2020; Var- 358 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards gas-Jaimes et al. 2021). The ensemble model (seven algorithms × 10 repetitions) was subsequently used to project the potential distributions of the species. To evaluate the potential distribution areas of Spiny Lizards, we applied two threshold criteria: the presence of training in the 10th percentile and the presence of training in the 80th percentile. The 10th percentile represents a more permissive threshold, often used to capture the full potential range of a species, including marginal habitats. In contrast, the 80th percentile is a more conservative and ecologically realistic threshold and is especially suitable for assessing vulnerability in the species studied. These analyses were performed using the EcoNicheS package in R. Landscape connectivity Functional connectivity analysis was performed with a circuit theory (CT) model, utilising a landscape scenario as a resistance surface (Mateo-Sánchez et al. 2014) implemented through Circuitscape v. 4.0 software (McRae et al. 2013). The resistance surface layers were derived from habitat suitability maps generated with the ENMeval package and computed using the inverse transformation equation (1 − MNE), where MNE represents the ecological niche model output. Given the constraints in modifying processor numbers within the Circuitscape v. 4.0 interface, Julia 1.8.2 (Bezanson et al. 2017) software was used to calculate habitat connectivity for species with abundant presence records. The resulting maps were visualised and exported in ArcMap. Furthermore, the identification of specific connectivity areas (corridors, linkages, stepping stones) in geographic space involved analysing the overlap between the Protected Natural Areas polygons and species-specific connectivity models. Habitat fragmentation analysis (number of patches) To quantify habitat fragmentation for each species, Sceloporus, we used threshold binary distribution maps with a presence value ≥ 80 (on a 0–100 suitability scale). Raster maps were processed using the terra package (Hijmans 2025) in R. Each binary map was first converted into a clump map, which assigns unique IDs to contiguous patches (clumps) of presence cells. The number of patches per species was then calculated using the clump output, excluding the NA values. This approach allowed us to estimate both the total number of habitat patches (a proxy for fragmentation) and the mean patch size by dividing the total area by the number of patches. All raster analyses were conducted at a consistent resolution and coordinate reference system across the study area. Anthropogenic pressure index and habitat exposure To quantify species exposure to anthropogenic pressures, we overlaid each species’ binary habitat suitability map with three rasterised land-use layers representing major human activities: urban areas (U), agriculture (A) and cultivated pastures (P). A threshold was applied for all habitat suitability maps using a value of ≥ 80 (on a 0–100 scale), producing binary maps of suitable habitat. These binary maps and land-use layers were processed in R using the terra package. Each anthropogenic pressure layer was coded as 1 (presence of pressure) or 0 (absence) and combined with the species’ distribution raster. 359 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards For each pixel within a species’ distribution, we calculated an Anthropogenic Pressure Index (API) using the following equation: API = U + A + P where: U = 1 if the pixel overlaps with urban areas, 0 otherwise; A = 1 if the pixel overlaps with agricultural land, 0 otherwise; P = 1 if the pixel overlaps with cultivated pastures, 0 otherwise. This pixel-level index ranges from 0 (no anthropogenic pressure) to 3 (overlap with all three pressures), representing the cumulative number of anthropogenic pressures affecting each cell. To summarie the overall impact, we calculated the percentage of suitable habitat under anthropogenic pressure for each species using the following formula: Pressure Index(%) = (Total suitable habitat cells) × 100 This index reflects the proportion of each species’ predicted habitat affected by at least one anthropogenic pressure. The results were visualised using horizontal bar plots with a continuous colour gradient (blue to red) to represent increasing levels of pressure, created using the ggplot2 package. This index was a key component of our integrative vulnerability assessment, helping to identify species whose ranges are disproportionately exposed to human landuse conversion across TMVB. Multivariate pressure zonation (PCA and k-means) To spatially delineate regions of cumulative anthropogenic pressure, we performed a principal component analysis (PCA) on the binary urban, agricultural and pasture layers using the RStoolbox (Müller et al. 2024) and stats (R Core Team 2024) packages. The first three principal components were extracted to represent key gradients in human land use. These components were then clustered using the kmeans function from base R (stats package) with k = 4 to define multivariate pressure zones across the TMVB. Each zone was assigned to one of four pressure classes: low pressure: natural or semi-natural areas. Moderate pressure: Low-intensity agriculture or pastures. High pressure: Mixed agricultural–urban mosaics. Very high pressure: Urban or industrialised areas. This zonation was used to calculate the distribution of each species across pressure categories. Composite vulnerability and pressure integration To assess overall extinction risk, we integrated two key metrics: the spatial vulnerability index (based on habitat fragmentation and isolation) and the total anthropogenic pressure index (based on the overlap with land-use threats). Both indices were plotted in a bivariate space to classify each species into one of four risk categories: High risk – high vulnerability and high exposure; Pressure-only – high exposure, but low spatial vulnerability; Vulnerability-only – high spatial fragmentation, but low exposure to anthropogenic pressure; and Low risk – low values for both indices. This classification was visualised using 366 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards the connectivity values are concentrated in the Jalisco Valley area and the central part of Nayarit State. Notably, a corridor spanning between these regions, particularly evident in the Municipality of Ixtlán del Río in Nayarit, facilitated connectivity for species such as S. albiventris, S. bulleri, S. horridus and S. utiformis. Habitat fragmentation, spatial vulnerability and composite vulnerability indices Binary distribution maps (threshold ≥ 80) revealed substantial variation in habitat size and fragmentation across species (Fig. 5). The number of habitat patches ranged from a single patch (e.g. S. scalaris) to more than 500 (e.g. S. minor) (Fig. 5A), indicating substantial heterogeneity in landscape continuity. Species with fewer and larger patches (e.g. S. melanogaster, S. mucronatus, S. scalaris and S. variabilis) exhibited lower fragmentation (Fig. 5A, B), while those with many small patches (e.g. S. minor, S. albiventris and S. anahuacus) were considered more fragmented (Fig. 5A, B). Five species (S. minor, S. albiventris, S. anahuacus, S. subniger and S. bulleri) exhibited the greatest vulnerability. The composite vulnerability index (Fig. 5C) integrates habitat area, number of habitat patches and average patch size. This finding revealed substantial variation amongst the Sceloporus species. Sceloporus minor exhibited the highest vulnerability score, indicating a combination of small, fragmented and isolated habitat patches. Furthermore, S. anahuacus, S. albiventris and S. megalepidurus also exhibited high vulnerability values, suggesting a concerning spatial configuration of their potential habitats that may restrict population viability and dispersal. Mid-vulnerability level scores were calculated for S. subniger, S. palaciosi and Table 3. Potential distribution area of lizards of the genus Sceloporus in the Mexican territory calculated in km2 for the two algorithms used, using the 10th percentile value. Species 10th percentile value ENMeval biomod2 km2km2 S. albiventris 0.166 112,705.6 187,709.32 S. anahuacus 0.167 11,994.5 57,995.51 S. bicanthalis 0.176 14,938 51,896.62 S. bulleri 0.384 97,380.2 82,466.86 S. dugesii 0.309 93,919.8 154,462.6 S. grammicus 0.211 162,105 280,138.91 S. horridus 0.239 163,002.2 217,905.88 S. megalepidurus 0.272 12,374.1 24,725.56 S. melanogaster 0.261 179,567 380,387.21 S. minor 0.235 127,956.5 270,620.19 S. mucronatus 0.175 73,933.7 457,427.91 S. palaciosi 0.289 14,326.2 14,359.08 S. scalaris 0.169 251,970.4 400,727.4 S. subniger 0.155 34,206.1 89,236.87 S. torquatus 0.189 28,902.4 52,763.1 S. utiformis 0.277 113,330.9 148,089.93 S. variabilis 0.198 267,366.8 150,5482.92 367 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards Figure 4. Landscape connectivity maps of 17 species of the genus Sceloporus. The gradient in the upper right corner indicates the degree of connectivity of each species in its study area. In the upper right corner, a map of Mexico (in grey) shows the study area M for each species, with the Trans-Mexican Volcanic Belt (TMVB) highlighted in green. S. bicanthalis, which may be at risk if current anthropogenic pressures persist or intensify. Conversely, species such as S. variabilis, S. scalaris, S. mucronatus and S. melanogaster had notably negative index values, indicating broader, more 368 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards continuous and less fragmented potential distributions. Our results highlight key differences in spatial vulnerability amongst endemic species in the TMVB and underscore the necessity for species-specific conservation strategies. Anthropogenic pressure indices by species Fig. 6A displays the total anthropogenic pressure indices (%) for each of the 17 analysed Sceloporus species in the Trans-Mexican Volcanic Belt (TMVB). Notably, Sceloporus megalepidurus and S. bicanthalis exhibited the highest levels of exposure (> 17%) to potential habitats overlapping with anthropogenic land uses (urban areas, cultivated pastures and agriculture). These species are closely related to the pressure indices of S. melanogaster, S. mucronatus and S. scalaris (~ 15%). In contrast, S. palaciosi, S. bulleri and S. minor exhibited the lowest levels of anthropogenic pressure, with values less than 10%, suggesting that their suitable habitats are currently less impacted by land conversion. Human-driven threats undermine species conservation through habitat degradation, fragmentation and loss of ecological connectivity – pressures that are magnified when ranges are small or when representation in protected areas is low (Bolom-Huet et al. 2022). Partitioning each species’ suitable habitat (threshold ≥ 80/100) into four land-use pressure classes revealed a consistent pattern: urban areas dominated the habitat composition for nearly Figure 5. Habitat fragmentation metrics by species. A. Number of habitat patches per species; B. Mean patch size (km2) across species; C. Composite vulnerability index based on normalised values of habitat area, patch number and patch size. The bars represent relative spatial vulnerability; red tones indicate higher risk due to fragmentation and isolation. 369 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards all Sceloporus species (Fig. 6B). For several species – S. megalepidurus, S. palaciosi, S. bicanthalis and S. anahuacus – more than 80% of the predicted suitable habitat overlapped with urban landscapes, indicating extremely high exposure to human modification; agricultural land generally comprised a smaller fraction and induced pastures contributed only marginally (Fig. 6B). For S. megalepidurus, S. palaciosi, S. bicanthalis and S. anahuacus, more than 80% of their habitat is located in urbanised areas, indicating extremely high exposure to human-modified landscapes (Fig. 6B). In contrast, agriculture (blue) contributed a smaller proportion of the total pressure, typically representing less than 20% of each species’ habitat (Fig. 6B). Only a few species, such as S. minor, S. horridus and S. dugesii, exhibit modest levels of agricultural overlap. Induced grasslands (Fig. 6, orange), which often emerge from land degradation or pasture conversion, occupy a small proportion of the habitat, primarily affecting S. palaciosi, S. megalepidurus and S. anahuacus. The absence of natural areas (Fig. 6, grey bars) for nearly all the species indicates that virtually the entire suitable habitat of these lizards is subjected to different levels of anthropogenic pressure (Fig. 6B). This aligns with earlier findings of high fragmentation and low protection coverage, supporting the urgency for conservation actions. The bivariate plot combines the vulnerability index and the total index of anthropogenic pressure. Fig. 6C shows specific conservation profiles amongst the 17 Sceloporus species analysed. Sceloporus megalepidurus is a “high risk” species because it is the only Figure 6. Anthropogenic pressure and vulnerability patterns across species. A. Total anthropogenic pressure index (%) for each species, derived from overlap with urban, agricultural and pasture layers; B. Proportional habitat exposure to four land-use pressure zones: agricultural (blue), urban (red), induced grassland (orange) and unpressured habitat (grey); C. Bivariate plot of vulnerability vs. total anthropogenic pressure. Species are colour-coded by risk category: red = high risk, orange = pressure only, yellow = vulnerability only and green = low risk. 370 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards species that occupies the upper-right quadrant, as indicated by high values for both vulnerability and pressure. Ten species – including S. mucronatus, S. melanogaster, S. scalaris and S. variabilis – are “pressure only” species because they are highly exposed to anthropogenic threats despite lower spatial vulnerability scores. In contrast, the lizards, S. minor and S. palaciosi, are “vulnerability only” species, characterised by fragmented or limited distributions with minimal overlap with human-modified landscapes. The remaining species, such as S. bulleri, S. torquatus and S. horridus, clustered in the lower-left quadrant, indicating low overall risk due to both low vulnerability and low exposure to anthropogenic pressures (Fig. 6C). Multivariate anthropogenic pressure analysis and pressure zonation (K-means classification) The analysis identified four distinct anthropogenic pressure zones across the Trans-Mexican Volcanic Belt (TMVB) (Fig. 7): Zone 1 (green) – Predominantly natural or semi-natural landscapes, including forests and high-elevation areas. These regions exhibit low anthropogenic pressure, but are sparsely distributed and highly fragmented, underscoring the limited extent of pristine habitats within the TMVB. Zone 2 (yellow) – Moderately pressured landscapes characterised by low-intensity land use, such as dispersed agriculture or pasture. These areas function as transitional buffers around urban zones and regions of intensive agriculture. Zone 3 (orange) – Mixed-use landscapes dominated by agricultural activity interspersed with urban development. This zone covers the majority of the TMVB, forming the primary habitat matrix in which most species must persist and disperse, highlighting the prevalence of human-modified environments in central Mexico. Zone 4 (white) – Areas of very high anthropogenic pressure, primarily urban and industrial zones concentrated around major metropolitan centers such as Mexico City, Toluca, Puebla and Guadalajara. These highly transformed environments offer minimal habitat suitability for native biodiversity. This zonation provides a spatial framework for assessing the cumulative impacts of land-use change (LUC) and identifying regions where conservation and restoration interventions are most urgent. The weighted vulnerability richness analysis revealed clear spatial patterns of conservation priority. The vulnerability hotspot map, generated by integrating binary species distributions (threshold ≥ 80) weighted by species-specific vulnerability indices, highlighted several high-priority regions. Hotspots – represented by brighter yellow and orange tones – were concentrated in the eastern TMVB (Puebla and Veracruz highlands), the Sierra de las Cruces between Mexico City and Toluca and montane areas of western Jalisco and Nayarit. Conversely, darker areas corresponded to regions with low species richness or a predominance of species with low vulnerability indices. This mapping approach offers powerful visualisation of conservation urgency, enabling the identification of priority areas where protection or restoration actions could simultaneously benefit multiple at-risk species. Temporal trend of forests cover loss (2001–2023) To assess deforestation patterns within the Trans-Mexican Volcanic Belt (TMVB), we analysed annual forest cover loss from 2001 to 2023 using data extracted from the Global Forest Watch platform (Fig. 8A). The dataset revealed 371 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards Figure 7. A. Multivariate pressure zonation across the TMVB. K-means clustering (k = 4) of PCA-derived axes from binary urban, agricultural and pasture layers. Green: low anthropogenic pressure (natural/semi-natural). Yellow: moderate pressure (low-density agriculture/pasture). Orange: high pressure (mixed agriculture-urban). White: very high pressure (urban/industrial). This zonation highlights the spatial structure of cumulative human impact; B. Vulnerability hotspot map (weighted richness). Spatial distribution of weighted species richness across the TMVB, with brighter areas indicating greater concentrations of vulnerable species. Vulnerability was weighted by each species’ composite index. Darker zones represent either low richness or species with low vulnerability scores. Figure 8. Temporal trend in forest loss across the TMVB (2001–2023). A. Annual forest loss (km2) from Global Forest Watch. The values indicate areas with > 30% canopy loss for trees > 5 m in length; B. Linear regression showing a statistically significant positive trend (R2 = 0.602, p = 1.35 × 10−5), indicating increasing deforestation over time. Point colours reflect loss intensity according to the viridis scale. Figure 9. Overlap between species distributions and protected natural areas. Bar chart showing the percentage of each species’ range that overlaps with federal (blue) and state-level (red) protected areas. The results highlight discrepancies in protection coverage amongst species, most of which rely heavily on state-level reserves. 372 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards fluctuations in tree cover loss across years, ranging from 12.2 km2 in 2002 to a peak of 127.0 km2 in 2018, with additional high-loss years observed in 2023 (94.2 km2) and 2019 (80.0 km2). The cumulative forest loss over the period was approximately 1242.0 km2 (Fig. 8A). These values correspond to areas with > 30% canopy loss for trees taller than 5 m and represent deforestation events potentially affecting Sceloporus habitats. A linear regression model was used to evaluate the temporal trend in forest cover loss (Fig. 8B). The results showed a statistically significant positive trend (R2 = 0.602, p = 1.35 × 10−5), indicating that forest loss increased considerably over time across the TMVB (Fig. 8B). These findings suggest that anthropogenic drivers of deforestation, such as urban expansion, agriculture, and fire, may be intensifying in this region. The trend was visualised with a scatter plot and a dashed regression line with 95% confidence intervals. Forest loss values were coloured using a viridis palette to enhance contrast between lowand high-loss years and each value was labelled directly on the graph. The pronounced upward trajectory reinforces the urgency of integrating deforestation dynamics into biodiversity vulnerability assessments and regional conservation strategies. Overlap with protected natural areas The species Sceloporus palaciosi, S. bicanthalis and S. anahuacus had the highest total protection, with up to ~ 33% of their distributions falling into protected natural areas and splitting between federal and state jurisdictions (Fig. 9). On the other hand, the species at the top of the graph (e.g. S. albiventris and S. utiformis) exhibited less than 20% protection (Fig. 9). In addition, the majority of their coverage is related to State-level protected areas, suggesting weak federal representation. A few species (e.g. S. melanogaster, S. scalaris and S. mucronatus) exhibited relatively balanced contributions from both protection types. Federal protection alone is rarely sufficient; most species rely on a combination of federal and state-level coverage to reach even minimal representation (Fig. 9). Discussion Understanding the environmental constraints and anthropogenic threats related to the distribution of reptile species is essential for developing conservation strategies. They must reflect the current and future pressures acting against biodiversity. In the case of the genus Sceloporus, a highly diverse and predominantly endemic group from Mexico, evaluating vulnerability beyond geographic range size has become increasingly important. Our study advances this goal by integrating multiple spatially explicit analyses considering habitat conditions and the intensity of human-induced landscape changes. A major strength of our research lies in the incorporation of land-use variables into the ecological niche modelling framework. While climatic variables (including temperature, precipitation and solar incidence) exhibited substantial contributions across most species, we also observed the significance of land-cover types in certain cases. Our findings align with those of Urbina-Cardona et al. (2006), who suggested a deep association between amphibian and reptile species and environmental factors such as vegetation cover, temperature and relative humidity, which indirectly influence factors, such as food and shelter availability. From the 373 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards models derived, two overarching distribution patterns within this species group can be delineated. Abies, species exhibiting broad potential distributions; positively correlated with urban areas (e.g. S. albiventris, S. dugesii, S. grammicus, S. horridus, S. mucronatus and S. variabilis); or with climatic variables (e.g. S. bulleri, S. melanogaster, S. minor, S. scalaris and S. utiformis), which extend across highland regions and the Eastern or Western Sierras. Conversely, species with more restricted distributions, primarily confined to the elevated areas of the TMVB (e.g. S. subniger, S. anahuacus, S. bicanthalis, S. megalepidurus, S. palaciosi and S. torquatus), are often associated with native vegetation, particularly Pinus and Abies forests. The negative relationship between the percentage of arid vegetation and the distributions of several species can be explained by their reliance on solar radiation to reach basal body temperatures. However, excessive exposure to high temperatures poses the risk of overheating and mortality, underscoring the importance of thermal shelters (Huey and Stevenson 1979; Sinervo et al. 2010). Notably, species with broader distribution ranges have demonstrated greater capacities for thermal acclimatisation than do those with more restricted ranges. This adaptive flexibility potentially enhances the ability of these Spiny Lizards to colonise new habitats and renders them less vulnerable to the adverse impacts of climate warming (Deutsch et al. 2008; Sunday et al. 2011; Markle and Kozak 2018). The preference for urban areas amongst certain species may be attributed to their adaptability to these environments, where they can find thermal shelters and other environmental amenities, albeit at the expense of associated drawbacks. Notably, the areas exhibiting the highest habitat suitability and connectivity values for several target species are situated within the mountainous regions of the TMVB, nested amidst sprawling metropolitan cities, such as Mexico City, the Toluca Valley and Puebla (See Suppl. material 1: fig. S2). The forested expanses of this region hold significant commercial value owing to their timber resources (FAO 2010; Pérez-Miranda et al. 2011) and the availability of water resources (Del Campo et al. 2019). Consequently, substantial deforestation events and land-use changes related to agricultural expansion and infrastructure development have been observed in recent decades (Velázquez et al. 2002; Rosete-Vergés et al. 2014; Monjardín-Armenta et al. 2017; CONAFOR 2020; De Alba Rosano et al. 2020), with indications that this trend will persist. The TMVB is a significant biodiversity hotspot in the country (Mastretta‐ Yanes et al. 2015; González-Fernández et al. 2018; Sunny et al. 2022; Rubio-Blanco et al. 2024). In this study, three pivotal sites for lizard distribution and connectivity were pinpointed. The pivotal sites are important for other endemic reptile species in the country (Sunny et al. 2023). The three sites included “Sierra de las Cruces”, which delineates the Toluca and Mexico City valleys and emerged as the primary connectivity hub for the seven species in our study. Despite the presence of federally protected natural areas such as “Desierto de los Leones” and “La Marquesa National Park”, along with State-level protected areas, such as “Otomí-Mexica Park”, these areas are experiencing significant degradation. Issues such as the expansion of upscale residential areas and associated infrastructure, similar to the “Toluca-Naucalpan highway”, are grave concerns (Mejía 2012). Strengthening connectivity between the “Nevado de Toluca”, the “Sierra de las Cruces” and the “Chichinautzin Biological Corridor” in the southern Toluca Valley region is deemed essential, given the moderate to 374 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards low connectivity levels indicated by our models. The valley area of Jalisco, bordering the State of Nayarit (See Suppl. material 1: fig. S3), emerged as the primary region for seven species, despite the presence of protected areas such as the “La Primavera” and “Sierra de Quilla” APFFs. In Nayarit, the Municipalities of Tepic, Santa María del Oro, Jala and Ixtlán del Río stand out due to their rugged relief, providing relevance for at least three species. Notably, a low-altitude relief pass was identified between “Sierra del Pinabete” and “Sierra Madre del Sur” in the Municipality of Ixtlán del Río, connecting it with the aforementioned area of Jalisco (INEGI 2005). Hence, this region is suggested to constitute a potential corridor for species connectivity in the western region. A comparison of our results with those of a previous study utilising the same algorithm (ENMeval) and environmental variables (Vargas-Jaimes et al. 2021) revealed that the calculated potential distribution areas for S. anahuacus, S. bicanthalis, S. megalepidurus and S. palaciosi were smaller than the area calculated for the endemic salamander Pseodoeurycea leprosa (34,824 km2). Despite being classified as “least concern” (LC) on the IUCN Red List and “Threatened” (A) according to NOM-059 and SEMARNAT (2002), due to its high association with forest cover variables and low temperatures, this salamander has a greater conservation status than these Spiny Lizards, suggesting a need for risk re-categorisation. The potential distributions of S. subniger and S. anahuacus are largely influenced by the presence of Abies forests, which are predominantly found in the mountainous regions of the TMVB and are increasingly threatened by encroaching on legal and illegal logging (González-Fernández et al. 2022). With rising global temperatures, shifts in species altitudinal ranges are anticipated (Bonino et al. 2014; Sillero 2021). Evidence of such shifts has already been observed in S. megalepidurus (Díaz de la Vega-Pérez et al. 2022) and S. utiformis (Muñoz-Nolasco et al. 2022), attributed to climate change-induced habitat alterations. Oviparous species of the genus Sceloporus are reportedly more vulnerable to global warming than sympatric viviparous species (Ma et al. 2022). However, Chiu-Valderrama et al. (2022) documented a contraction in the distribution range of S. grammicus (a viviparous species) under climate change scenarios. To mitigate the effects of isolation and limited connectivity, we propose the implementation of local management strategies. These strategies involve enriching the environment by increasing thermal and ecological refuges in areas devoid of vegetation, such as isolated trees located at the intersection of forest fragments, conserving secondary vegetation along forest edges and establishing “living fences” (Urbina-Cardona et al. 2006). Additionally, measures should be taken to avoid the extraction of stones, which serve as important shelter and basking sites. There is also an urgent need to address deforestation practices and prevent the invasion of exotic domestic species in natural environments (Sunny et al. 2024b). For species associated with anthropogenic environments, strategies should focus on fostering a harmonious reptile–human relationship and dispelling misconceptions and aversive beliefs regarding this group of organisms. A key objective of our study was to compare Maxent models optimised with ENMeval with ensemble models optimised with biomod2. Various studies have demonstrated that models with different levels of complexity can yield significantly different predictions of suitable areas for species (Elith et al. 2010; Anderson and González 2011; Shcheglovitova and 375 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards Anderson 2013). Although the evaluation values of the models produced by both algorithms indicated high reliability in terms of the predictions, the Ensemble models proposed broader distribution ranges than those generated by ENMeval for more than half of the species, as indicated by the calculated area. To increase the accuracy of ecological niche estimations for these and other species, it is essential to incorporate variables that capture interactions within the biotic space, such as those representing co-habiting species in the geographical area (Araújo and Luoto 2007). Furthermore, variables indicating human impact should be included in estimating suitability gradients and tolerance to disturbance (Ochoa-Ochoa and Ríos-Muñoz 2019). The potential distribution models proposed in this study represent valuable tools for conservation efforts and decision-making processes. Regarding the ecological niche modelling, we observed a significant contribution of climatic variables, particularly maximum temperature and solar radiation level, which aligns with expectations given the ectothermic nature of these organisms (Rivera-Rea et al. 2023; Plasman et al. 2025). However, certain land-cover variables, notably Pinus forests, Abies forests and urban areas, also exhibited notable contributions to the presence of several species. TMVB has emerged as a critical area for the distribution and conservation of Spiny Lizards, with projections indicating high habitat suitability and connectivity for most species within this region. Within the TMVB, we identified three sites characterised by these favourable conditions, with the Sierra de las Cruces standing out as particularly crucial for the conservation and maintenance of connectivity amongst several species studied here, as well as other endemic species in the area. Additionally, we recommend revisiting the risk categorisation of certain species, such as S. anahuacus, S. bicanthalis, S. megalepidurus and S. palaciosi. The spatial prioritisation analysis integrating vulnerability and exposure to anthropogenic pressures revealed that five species – S. anahuacus, S. albiventris, S. bicanthalis, S. megalepidurus and S. minor – exhibited both high habitat fragmentation and low protection coverage, classifying them as high priorities for conservation. Our analysis of habitat composition across multivariate pressure zones revealed an alarming pattern: the vast majority of suitable habitats for Sceloporus species are embedded within anthropogenically modified landscapes, particularly urban areas. Urban zones account for more than 80% of the predicted habitats of several species, including S. megalepidurus, S. palaciosi and S. bicanthalis, all of which are endemic to highland regions of the TMVB. These results suggest a substantial mismatch between current land-use dynamics and the ecological requirements of these lizards. The nearly complete absence of habitat in “no pressure” zones for most species underscores a critical conservation concern: even species currently classified as “least concern” (e.g. S. anahuacus and S. grammicus) may face hidden risks due to high exposure to urbanisation and habitat transformation. The composite vulnerability index also provided key insights. Species with the highest vulnerability scores, S. minor, S. anahuacus, S. albiventris and S. megalepidurus, were characterised by highly fragmented and reduced habitat patches, posing threats to meta-population dynamics and long-term viability. Conversely, species such as A, S. scalaris and S. mucronatus had low vulnerability scores, likely reflecting broader, more continuous habitats that buffer them against anthropogenic disturbance. Bivariate 382 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards Spain). International Journal of Geographical Information Science 28(8): 1531–1546. https://doi.org/10.1080/13658816.2014.889299 McRae BH, Shah VB, Mohapatra TK, Anantharaman R (2013) Circuitscape 4. Seattle WA: The Nature Conservancy. Mejía HB (2012) Deterioro socioambiental en la Sierra de las Cruces. Diseño y Sociedad, (32): 40–49. Metz CE (1978) Basic principles of ROC analysis. Seminars in Nuclear Medicine 8(4): 283–298. [WB Saunders] https://doi.org/10.1016/S0001-2998(78)80014-2 Monjardín-Armenta SA, Pacheco-Angulo CE, Plata-Rocha W, Corrales-Barraza G (2017) La deforestación y sus factores causales en el estado de Sinaloa, México. Madera y Bosques 23(1): 7–22. https://doi.org/10.21829/myb.2017.2311482 Morrone JJ (2014) Biogeographical regionalisation of the Neotropical region. Zootaxa 3782(1): 1–110. https://doi.org/10.11646/zootaxa.3782.1.1 Müller K, Schwalb-Willmann J, Wegmann M, Leutner B (2024) RStoolbox: An R package for remote sensing data analysis. Methods in Ecology and Evolution 15(2): 345–358. https://doi.org/10.1111/2041-210X.14451 Muñoz-Nolasco FJ, Cruz-Sáenz D, Méndez-de la Cruz FR (2022) Extensión de distribución altitudinal superior de Sceloporus utiformis Cope, 1864 (Squamata: Phrynosomatidae) en la reserva de la biósfera sierra de Manantlán, México. Revista Latinoamericana de Herpetología 5(4): 1–5. https://doi.org/10.22201/fc.25942158e.2022.4.511 Muscarella R, Galante PJ, Soley-Guardia M, Boria RA, Kass JM, Uriarte M, Anderson RP (2014) ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods in Ecology and Evolution 5: 1198–1205. https://doi.org/10.1111/2041210X.12261 Ochoa-Ochoa LM, Flores-Villela O (2006) Áreas de diversidad y endemismo de la herpetofauna Mexicana. UNAM-CONABIO, México, D.F., 211 pp. Ochoa-Ochoa LM, Ríos-Muñoz CA (2019) Estimación de adecuabilidad de especies en ambientes antropizados desde el enfoque de modelado de nicho. Antropización: primer análisis integral (Ornelas-García, CP, F. Álvarez, and AL Wegier, Eds). Red Temática de Fauna Nativa en Ambientes Antropizados-Consejo Nacional de Ciencia y Tecnología, Instituto de Biología-Universidad Nacional Autónoma de México. Ciudad de México, México, 345–362. https://doi.org/10.22201/ ib.9786073020305e.2019.c18 Osorio-Olvera L, Lira-Noriega A, Soberón J, Peterson AT, Falconi M, Contreras-Díaz RG, Martínez-Meyer E, Barve V, Barve N (2020) ntbox: An R package with graphical user interface for modeling and evaluating multidimensional ecological niches. Methods in Ecology and Evolution 11(10): 1199–1206. https://doi. org/10.1111/2041-210X.13452 Pérez-Miranda R, Valdez-Lazalde JR, Moreno-Sánchez F, González-Hernández A, Valdez-Hernández JI (2011) Predicción espacial de cambios del uso de suelo en Texcoco, Estado de México. Revista Mexicana de Ciencias Forestales 2(5): 59–72. https://doi.org/10.29298/rmcf.v2i5.586 Peterson AT, Papeş M, Soberón J (2008) Rethinking receiver operating characteristic analysis applications in ecological niche modeling. Ecological Modelling 213(1): 63–72. https://doi.org/10.1016/j.ecolmodel.2007.11.008 Phillips SJ, Anderson RP, Schapire RE (2006) Maximum entropy modeling of species geographic distributions. Ecological Modelling 190(3–4): 231–259. https://doi. org/10.1016/j.ecolmodel.2005.03.026 383 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards Plasman M, Díaz de la Vega-Pérez AH, McCue MD, Tepos-Ramírez M, Reynoso VH (2025) The ultimate challenge to climate change: Endurance of a thermophilic reptile to the harsh temperatures on an extremely hot island. PLOS ONE 20(4): e0320796. https:// doi.org/10.1371/journal.pone.0320796 R Core Team (2024) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/ R Studio Team (2021) RStudio: Integrated Development Environment for R. RStudio, PBC, Boston, MA. Rivera-Rea J, Macotela L, Moreno-Rueda G, Suárez-Varón G, Bastiaans E, Quintana E, González-Morales JC (2023) Thermoregulatory behavior varies with altitude and season in the sceloporine mesquite lizard. Journal of Thermal Biology 114: 103539. https://doi.org/10.1016/j.jtherbio.2023.103539 Rosas-Espinoza VC, Rodríguez-Zaragoza FA, Álvarez-Grzybowska E, Peña-Joya KE, Santiago-Pérez AL, Godoy-González AA, Huerta-Martínez FM (2024) Taxonomic and functional diversity of reptiles in a heterogeneous landscape of Jalisco State, West-Central Mexico. Diversity 16(7): 394. https://doi.org/10.3390/d16070394 Rosete-Vergés FA, Pérez-Damián JL, Villalobos-Delgado M, Navarro-Salas EN, Salinas-Chávez E, Remond-Noa R (2014) El avance de la deforestación en México 1976–2007. Madera y Bosques 20(1): 21–35. https://doi.org/10.21829/ myb.2014.201173 Rubio-Blanco T, Martínez-Díaz-González R, Heredia-Bobadilla RL, Guido-Patiño JC, Arenas S, Caballero-Viñas C, Manjarrez J, Domínguez-Vega H, Ortiz YG, Ramos-Olguin AD, Sunny A (2024) Predicting the effects of climate and land use changes on small rattlesnakes in central Mexico: Insights for conservation planning. Journal for Nature Conservation 79: 126607. https://doi.org/10.1016/j.jnc.2024.126607 Ruiz-Reyes J, Heredia-Bobadilla RL, Ávila-Akerberg V, Tejocote-Perez M, Gómez-Ortiz Y, Domínguez-Vega H, Ramírez-Corona F, Alvarez-Lopeztello J, Bolom-Huet R, Sunny A (2024) Assessing functional connectivity and anthropogenic impacts on Ambystoma altamirani populations in Bosque De Agua, Central Mexico. European Journal of Wildlife Research 70(5): 85. https://doi.org/10.1007/s10344-024-01838-8 Segarra J, González-Fernández A, Osorno-Covarrubias J, Couturier S (2024) The role of critical remote sensing in environmental justice struggles. Progress in Environmental Geography 3(3): 185–211. https://doi.org/10.1177/27539687241269331 SEMARNAT (2002) [March 06] Norma Oficial Mexicana NOM-059-ECOL-2001, Protección ambiental – Especies nativas de México de flora y fauna silvestres – Categorías de riesgo y especificaciones para su inclusión, exclusión o cambio – Lista de especies en riesgo. Diario Oficial de la Federación. SEMARNAT (2025) Mexican Official Norm NOM-059-SEMARNAT-2025, Environmental protection – Native species of Mexico – Wildlife – Risk categories and specifications for inclusion, exclusion or change – and list of species at risk. Diario Oficial de la Federación, Mexico City, Mexico. Shcheglovitova M, Anderson RP (2013) Estimating optimal complexity for ecological niche models: A jackknife approach for species with small sample sizes. Ecological Modelling 269: 9–17. https://doi.org/10.1016/j.ecolmodel.2013.08.011 Sillero N (2021) Climate change in action: Local elevational shifts on Iberian amphibians and reptiles. Regional Environmental Change 21: 101. https://doi.org/10.1007/ s10113-021-01831-w Sinervo B, Méndez-de-la-Cruz F, Miles DB, Heulin B, Bastiaans E, Cruz MV-S, Lara-Resendiz R, Calderón-Espinosa NM-MML, Meza-Lázaro RN, Gadsden H, 384 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards Avila LJ, Morando M, Riva IJDl, Sepulveda PV, Rocha CFD, Ibargüengoytía N, Puntriano CA, Massot M, Lepetz V, Oksanen TA, Chapple DG, Bauer AM, Branch WR, Clobert J, Sites Jr JW (2010) Erosion of lizard diversity by climate change and altered thermal niches. Science 328(5980): 894–899. https://doi.org/10.1126/ science.1184695 Sites JW, Archie JW, Cole CJ, Flores-Villela O (1992) A review of phylogenetic hypotheses for Spiny Lizards of the genus Sceloporus (Phrynosomatidae): Implications for ecological and evolutionary studies. Bulletin of the American Museum of Natural History 213: 5–28. Smith GR, Ballinger RE (2001) The ecological consequences of habitat and microhabitat use in Spiny Lizards: a review. Contemporary Herpetology 1: 1–28. https://doi. org/10.17161/ch.vi1.11957 Smith H (1941) Las provincias bióticas de México, According to la distribución geográfica de las lagartijas del género Sceloporus. Anales. Escuela Nacional de Ciencias Biológicas (México) 2: 103–110. Sokal RR (1995) The principles and practice of statistics in biological research. Biometry, 451–554. Suazo-Ortuño I, Ramírez-Bautista A, Alvarado-Díaz J (2023) Amphibians and reptiles of Mexico: Diversity and conservation. In: Horváth C, López-Pujol J, Newmark WD (Eds) Mexican Fauna in the Anthropocene. Springer International Publishing, Cham, 105–127. https://doi.org/10.1007/978-3-031-17277-9_6 Sunday JM, Bates AE, Dulvy NK (2011) Global analysis of thermal tolerance and latitude in ectotherms. Proceedings of the Royal Society B 278(1713): 1823–1830. https:// doi.org/10.1098/rspb.2010.1295 Sunny A, Duarte-deJesus L, Aguilera-Hernández A, Ramírez-Corona F, Suárez-Atilano M, Percino-Daniel R, Manjarrez J, Monroy-Vilchis O, González-Fernández A (2019) Genetic diversity and demography of the critically endangered Roberts’ false brook salamander (Pseudoeurycea robertsi) in Central Mexico. Genetica 147: 149–164. https://doi.org/10.1007/s10709-019-00058-2 Sunny A, González-Fernández A, D’Addario M (2017) Potential distribution of the endemic imbricate alligator lizard (Barisia imbricata imbricata) in highlands of central Mexico. Amphibia-Reptilia 38(2): 225–231. https://doi.org/10.1163/15685381-00003092 Sunny A, López‐Sánchez M, Ramírez‐Corona F, Suárez‐Atilano M, González‐Fernández A (2022) Genetic diversity and functional connectivity of a critically endangered salamander. Biotropica 54(1): 42–56. https://doi.org/10.1111/btp.13025 Sunny A, Manjarrez J, Caballero-Viñas C, Bolom-Huet R, Gómez-Ortiz Y, Domínguez-Vega H, Heredia-Bobadilla RL, Torres-Romero EJ, González-Fernández A (2023) Modeling the effects of climate and land-cover changes on the potential distribution and landscape connectivity of three earth snakes (Genus Conopsis, Günther 1858) in central Mexico. Die Naturwissenschaften 110(6): 52. https://doi. org/10.1007/s00114-023-01880-7 Sunny A, Martínez-Valerio JC, Bolom-Huet R, Heredia-Bobadilla RL, Alvarez-Lopeztello J, Gómez-Ortiz Y, Guido-Patiño JC, Manjarrez J, González-Pedroza MG, Domínguez-Vega H (2024a) Modeling the present and future distribution of Ambystoma altamirani in the Transmexican Volcanic Belt, Mexico. Nature Conservation 56: 275–293. https://doi.org/10.3897/natureconservation.56.139402 Sunny A, Ruiz-Reyes J, Domínguez-Vega H, Gómez-Ortiz Y, Heredia-Bobadilla RL, Avila-Akerberg V, Manjarrez J, Reyes-Olivares E, García-Rendon S (2024b) Niche overlap by invasion of Oncorhynchus mykiss on the habitat of its amphibian prey 385 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards in central Mexico. Biological Invasions 26: 2183–2201. https://doi.org/10.1007/ s10530-024-03304-7 Sunny A, Marmolejo C, Vidal-López R, Falconi-Briones FA, Cuervo-Robayo ÁP, Bolom-Huet R (2025) EcoNicheS: Enhancing ecological niche modeling, niche overlap and connectivity analysis using the shiny dashboard and R package. PeerJ 13: e19136. https://doi.org/10.7717/peerj.19136 Thuiller W, Georges D, Gueguen M, Engler R, Breiner F, Lafourcade B, Patin R, Blancheteau H (2025) biomod2: Ensemble Platform for Species Distribution Modeling. R package version 4.2-6-2. https://CRAN.R-project.org/package=biomod2 Tischendorf L, Fahrig L (2000) On the usage and measurement of landscape connectivity. Oikos 90(1): 7–19. https://doi.org/10.1034/j.1600-0706.2000.900102.x Uetz P, Freed P, Aguilar R, Reyes F, Kudera J, Hošek J [Eds] (2025) The Reptile Database. http://www.reptile-database.org [Accessed on 09 November 2025] Urbina-Cardona JN, Olivares-Pérez M, Reynoso VH (2006) Herpetofauna diversity and microenvironment correlates across a pasture–edge–interior ecotone in tropical rainforest fragments in the Los Tuxtlas Biosphere Reserve of Veracruz, Mexico. Biological Conservation 132(1): 61–75. https://doi.org/10.1016/j.biocon.2006.03.014 Vargas-Jaimes J, González-Fernández A, Torres-Romero EJ, Bolom-Huet R, Manjarrez J, Gopar-Merino F, Pacheco XP, Garrido-Garduno T, Chávez C, Sunny A (2021) Impact of climate and land cover changes on the potential distribution of four endemic salamanders in Mexico. Journal for Nature Conservation 64: 126066 https://doi. org/10.1016/j.jnc.2021.126066. Velázquez A, Mas JF, Díaz JR, Mayorga-Saucedo R, Alcántara PC, Castro R, Fernández T, Bocco G, Escurra E, Palacio JL (2002) Patrones y tasas de cambio de uso del suelo en México. Gaceta Ecológica 62: 21–37. Venegas-Barrera CS, Manjarrez J, Rodríguez-Moreno Á, Mendoza-Walle YA, Horta-Vega JV, Rodríguez-deLeón IR, Sunny A, Azuara Domínguez A (2024) Representativeness, complementarity, and degree of local extirpation risk for Thamnophis species inside and outside of protected areas of Mexico. Ecologies 5(4): 697–715. https://doi. org/10.3390/ecologies5040041 Warren DL, Glor RE, Turelli M (2010) ENMTools: A toolbox for comparative studies of environmental niche models. Ecography 33: 607–611. https://doi.org/10.1111/ j.1600-0587.2009.06142.x Warren DL, Seifert SN (2011) Ecological niche modeling in Maxent: The importance of model complexity and the performance of model selection criteria. Ecological Applications : A Publication of the Ecological Society of America 21: 335–342. Wickham H (2016) ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag, New York. https://doi.org/10.1890/10-1171.1 Wiens JJ, Kuczynski CA, Arif S, Reeder TW (2010) Phylogenetic relationships of phrynosomatid Spiny Lizards based on nuclear and mitochondrial data, and a revised phylogeny for Sceloporus. Molecular Phylogenetics and Evolution 54(1): 150–161. https://doi.org/10.1016/j.ympev.2009.09.008 Wiens JJ, Reeder TW (1997) Phylogeny of the Spiny Lizards (Sceloporus) based on molecular and morphological evidence. Herpetological Monograph 11: 1–101. https://doi.org/10.2307/1467007 Zuur AF, Ieno EN, Walker NJ, Saveliev AA, Smith GM (2009) Mixed Effects Models and Extensions in Ecology with R. Springer, New York, 574 pp. https://doi.org/10.1007/9780-387-87458-6 386 Nature Conservation 59: 351–386 (2025), DOI: 10.3897/natureconservation.59.159660 Jossye Vargas-Jaimes et al.: Vulnerability assessment of Spiny Lizards Supplementary material 1 Supporting information of assessment of the vulnerability of Spiny Lizards Authors: Jossye Vargas-Jaimes, Javier Manjarrez, Fabiola Ramírez-Corona, Guillermo A. Woolrich-Piña, María G. González-Pedroza, María Elena Estrada-Zúñiga, Carlos Alejandro Rangel-Patiño, Crystian Sadiel Venegas-Barrera, Aurelio Nieto-Trujillo, Alicia Monserrat Vazquez-Marquez, Juan Carlos Guido-Patiño, Armando Sunnya Data type: docx Explanation note: Supporting information of Assessment of the vulnerability of Spiny Lizards in Central Mexico under anthropogenic pressure: integrating habitat suitability, landscape connectivity, fragmentation, and protected area coverage. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/natureconservation.59.159660.suppl1