Colombia is crossed by the Andes mountain range, which branches into three mountain ranges, giving it a high orographic complexity (Rangel, 2010). The Baudó system, Sierra Nevada de Santa Marta, and Macarena also stand out, complemented by valleys and mountains that confer a variety of meso- and microclimates (IDEAM-UNAL, 2018). Seventy percent of Colombia has a mean temperature of 24 °C (IDEAM-UNAL, 2018). The rainiest areas (Pacific region, Amazonian foothills, and plains) can have rainfall >4000 m/year, whereas the driest regions (Guajira Peninsula, some regions of the inter-Andean valleys) have rainfall of 500 m/year (IDEAM-UNAL, 2018). Six natural regions are recognized in Colombia: the Andean, Pacific, Amazon, Orinoco, Caribbean, and Insular regions (Fig. S1; Rangel, 2010).

Protected areas in Colombia began in 1938, consolidated with the signing of the CBD in 1992, ratified in Law 165 of 1994 (Lenis, 2014), and in 1993 the SINAP - National System of Protected Areas was proposed (Lenis, 2014). In 2019 in Colombia, OECMs were identified, strengthened, and reported at the international level (Santamaría et al., 2021). By 2023, Colombia had 1652 PA and 48 OECMs, encompassing 187,817 and 128,480 km2, respectively (SINAP, 2023; UNEP-WCMC, 2024; Table S1). The Protected Planet database (https://www.protectedplanet.net) was used as source of the PA and OECMs.

We calculated a remaining habitat map to evaluate the extent to which each species range was affected by habitat loss in 2022 (Appendix S1 in SM). Values of remaining habitat map closer to 1 indicate landscapes with a higher proportion of habitat (Table S2). We graphically explored the relationship between potential snake species richness (i.e., based on stacked species distribution models of semi-binary models; see details below) and remaining habitat both in geographic space and at the cell level (Velazco et al., 2023).

We obtained a list of 334 snake species native to Colombia from the Reptile DataBase portal (Uetz et al., 2023). We then reviewed the literature and other data sources for Colombia (BioModelos, 2023; SIB Col, 2023) to eliminate species with insufficient information on their occurrence. Species records were compiled from various open-access databases (Table S3). To better represent the environmental requirements of the species, models were constructed with records from throughout species' natural range (i.e., inside and outside Colombia). Taxonomic errors, geographic inaccuracies, and biases in record data lead to decreased performance of species distribution models (SDMs) and alterations in diversity patterns (Baker et al., 2022; Maldonado et al., 2015). We used the R bdc package (Ribeiro et al., 2022) to integrate record databases and perform spatial and temporal corrections (Appendix S2). Final database comprised 56,830 records and 309 species. Furthermore, for those species with no occurrences, <3 occurrences, or low model performance (i.e., Sorensen values < 0.7), but with distribution polygons available in Roll et al. (2017), we randomly sampled 1000 presences throughout species polygon.

SDMs were constructed using climatic and elevation environmental variables; however, some species were modeled using edaphic or hydrological variables, depending on their biological characteristics (i.e., aquatic, fossorial; Table S4). Thus, we initially chose 12 bioclimatic variables obtained from CHELSA v2.1 (Karger et al., 2017) and elevation (Jarvis et al., 2008), both with 1 km resolution. Edaphic variables were sourced from SoildGrid v.2.0 (Hengl et al., 2017) at 0−5 cm depth and 250 m resolution. We used Compound topographic index as a hydrological variable at 1 km resolution, as it serves as a proxy variable for watercourses (Table S4). The compound topographic index was obtained from the geomorpho90m database (Amatulli et al., 2020). All variables were upscaled to 5 km resolution with a geographic extent from the northern United States to southern South America. Initially, 17 variables were considered (Table S4). To reduce multicollinearity and the number of predictor variables, we constructed a Pearson correlation matrix (Fig. S2) using the values of the variables represented in all pixels. For all pairs of variables with a correlation ≥ |0.7|, we chose the variable with the highest biological significance. Finally, eight variables were selected (Table S4). Each species was fitted with a specific combination of variables (climatic, edaphic, hydrological, and elevation, Table S5). Although the selected variables were not correlated, the correlation structure could change for different species training areas (see species distribution models) (De Marco and Nóbrega, 2018). To overcome potential multicollinearity problems and reduce the number of predictors, we performed principal component analysis specific to each species variable combination (i.e., Table S4) and training area. Principal component analyses were performed based on correlation matrixes and selected a number of principal components that explained up to 95% of original variance (De Marco and Nóbrega, 2018). The derived principal components were used as predictors in SDMs.

We used SDMs to estimate species distribution and habitat suitability. SDMs predict species habitat suitability and geographic distribution by relating georeferenced observations (i.e., records) with environmental predictors (Soberón et al., 2017). The SDMs were created using eight algorithms (Appendix S3). We used several algorithms because no single algorithm can deal with all modeling conditions (e.g., species prevalence, niche breadth, or records number), and allows the use of consensus models (Qiao et al., 2015). We used flexsdm v1.3.6 R package to create SDMs, which allows the creation of flexible modeling protocols, and structuring functions in pre-modeling, modeling, and post-modeling steps (Velazco et al., 2022b).

Training areas of SDMs can affect environmental quality patterns and performance metrics (Barve et al., 2011). We delimited the training area of each species using a minimum convex polygon based on species occurrences plus a buffer of 500 km. We used a technique to filter records in environmental space to reduce the sampling bias of species with > 50 records (Appendix S3). Because we did not have absence data, we randomly sampled pseudo-absences distant to 50 km from records within the training area (Appendix S3). For species with occurrences sampled from their distribution polygons, we sampled pseudo-absences outside the species polygons (Mancini et al., 2024).

It is advisable to establish different modeling protocols depending on the requirements of the species and the number of records, as these affect the distribution pattern of the species and spatial prioritization analyses for conservation (Pimenta et al., 2022). Therefore, we designed three modeling protocols established based on records number (Appendix S4; Tables S5, S6, and S7).

We used Sorensen, Area Under the Curve (AUC), and True Skill Statistic (TSS) as model performance metrics (Fig. S3). The threshold that maximizes Sorensen’s metric was used to binarize models. The final model consisted of a consensus model calculated based on the arithmetic mean of environmental suitability values. For this procedure, we used only algorithms with Sorensen values ≥ 0.7.

SDM projections over large areas can lead to overprediction of high suitability zones outside the species’ current range, impacting spatial prioritization and diversity metrics (Velazco et al., 2020). To address this, we applied a minimum convex polygon with a 100 km buffer around it (Mendes et al., 2020) and used the threshold maximizing the Sorensen metric for binary species distributions.

Gap analysis assesses how well conservation areas meet representation targets for biological diversity (Rodrigues et al., 2004). Representation targets are based on species range size; species with ≤ 1000 km² require 100% of their range to be protected (i.e., within PA or OECMs), while those with ≥ 250,000 km² require 10%. Targets for species with ranges between 1000–250,000 km² were interpolated (Rodrigues et al., 2004). Species were classified as (1) 'Not Protected' if entirely unprotected, (2) 'Gap' if < 20% of the target is met, (3) 'Partial GAP' if 20-90% is met, and (4) 'Protected' if > 90% is met (Frederico et al., 2018; Table S8).

To identify priority areas for snake conservation, we used Zonation v4 (Moilanen et al., 2014), which ranks landscape cells based on their importance for conservation using various cell removal rules, weightings, and constraints (Di Minin et al., 2014). Zonation generates a top-down ranking of all cells in the study area, based on complementarity and irreplaceability (Moilanen et al., 2005). We selected two cell removal rules, core-area Zonation (CAZ), which prioritizes cells with rare or highly weighted species, and Additive Benefit Function (ABF); which emphasizes species richness (Di Minin et al., 2014); thus, both rules provided complementary results. Species were weighted according to their degree of endemism and threat categories (Appendix S5).

We used semi-binary models (i.e., environmental suitability values greater than the threshold were kept continuous while lower ones were set to zero) to reduce the problems of inflating the spatial prioritization analysis with many cells with low environmental suitability (Domisch et al., 2019). We used habitat loss as a cost layer, calculated as the inverse of remaining habitat (Appendix S6). Prioritization solutions were categorized into four classes by selecting 5%, 10%, 20%, and 30% of the highest priority cells and calculated proportion of priority cells within PA and PA + OECMs. The selection of 30% was motivated by the 30 × 30 target (CBD, 2021).

Colombia is recognized as a megadiverse territory and one of the 14 countries with the highest biodiversity worldwide (Correa, 2011). The Pacific and Andean regions had the highest snake diversity, which also stand out for their high diversity of other vertebrate groups (e.g., birds and amphibians) (Pinto-Erazo et al., 2020; Vélez et al., 2021). Although biogeographic processes explain the overall diversity of the Neotropics and the Andes (Turchetto-Zolet et al., 2013), snake richness patterns found in Colombia are understudied. We infer that snake richness in the Andean region may be driven by topographic and climatic heterogeneity that characterizes this region (Rangel, 2015, Rangel, 2010), as many studies have found a positive relationship between environmental heterogeneity and species richness (e.g., Körner, 2000; Stein et al., 2014). In contrast, high species richness in the Pacific region may be related to climatic stability (Rangel and Arellano, 2004). Pacific region has a low seasonality of precipitation and temperature, is characterized as one of the rainiest places worldwide (Rangel and Arellano, 2004) and harbors high species endemism (Rangel, 2015). Similar to other groups, the Amazon region also has considerable snake richness (Jenkins et al., 2013; Zapata-Ríos et al., 2022). It is presumed that the environmental heterogeneity associated with the formation of the Andes and the fluctuation of seasonal flooding in large alluvial river floodplains could have promoted Amazon diversity (Zapata-Ríos et al., 2022). Orinoquía and Caribbean regions exhibit lower snake richness, despite their high plant and vertebrate diversity (Jenkins et al., 2013). The lower richness towards the eastern region of Colombia appears to be consistent with other reptile and amphibian groups (IUCN, 2022a,b; Roll et al., 2017). However, the Colombian Amazon and Orinoquía present scarce species records (Suárez et al., 2021). Such data scarcity may be due to difficult access and lack of road infrastructure, in addition to being strongly affected by more than five decades of armed conflict (Cairo et al., 2018). Further studies are necessary to clarify the environmental and historical factors influencing these diversity patterns.

Colombia's PA and OECMs do not align with the highest snake richness areas, which are mostly outside these conservation networks (Figs. 1, 2). Such mismatch is consistent with amphibians and reptiles in some Colombian regions (Calderón-Caro et al., 2022; Valencia-Zuleta et al., 2014). Generally, PA suffer from different spatial bias types (Baldi et al., 2017), leading to low PA connectivity and poor representation of different ecosystems (Saura et al., 2017). Typically, PA are created in areas with touristic value, unproductive, or remote locations (Baldi et al., 2017; Phillips, 2007). Particularly Colombia's PA network started with the consolidation of some areas of economic interest, such as the Valle del Cauca, to support the sugarcane industry (Lenis, 2014). Although countries such as Colombia and Indonesia have described the limitations and challenges involved in OECMs (Atehortúa Arredondo et al., 2023; Estradivari et al., 2022), there is currently no information on their possible biases or effectiveness. However, expert consensus shed light on the opportunities and challenges of implementing OECMs (Alves-Pinto et al., 2021; Maini et al., 2023).

It has been repeatedly found that PA alone are not effective in representing biodiversity, as seen globally in mammals (Brum et al., 2017) and herpetofauna (Sánchez-Fernández and Abellán, 2015). In the case of Colombian snakes, it was no different. However, OECMs with PA doubled the representation of snake species richness. The GAP analysis showed that the categories "Partially GAP" and "Protected" also increased their representation target by 20.3% and 5.4%, respectively, for PA + OECMs. Therefore, OECMs play a crucial role in achieving conservation goals in Colombia (Rodríguez-Rodríguez et al., 2021). At the global level, OECMs can increase habitat connectivity and biodiversity representativeness (Alves-Pinto et al., 2021).

In South America, agricultural areas have expanded rapidly over the last 30 years (Eva et al., 2004; Sy et al., 2015). For example, Brazil lost the most forest cover from 1985 to 2004, mainly due to agriculture (Sy et al., 2015; Zalles et al., 2021). Bolivia and Peru have increased the areas of anthropic uses (agriculture, livestock, and urbanization) (MapBiomas, 2023, 2024). Although Colombia created policies focused on the management and sustainable land-use (MinAmbiente, 2016, 2013), it is not different from other countries in the region, as it has lost 7.4% of its natural cover in the last three decades (Fundación Gaia Amazonas, 2023).

When exploring the patterns of the remaining habitat with the patterns of species richness, we found that the areas with the highest species richness were consistent with the areas with the lowest remaining habitat, particularly in the Caribbean and parts of the Andean region (Fig. 4a). Frequently, regions with low topographic heterogeneity, flat, coastal, and low-elevation areas have the greatest urban and agricultural development (Gao et al., 2023; Rose et al., 2023). The Caribbean region is mainly composed of low and flat lands, and ports, industrial, and agricultural activities are developed in this region (Meisel-Roca and Pérez-Valbuena, 2006). The inter-Andean valleys of the Andean region are mainly used for agriculture (DANE, 2024). Although we found no clear pattern between species richness and remaining habitat, some studies have highlighted a positive relationship (Suazo-Ortuño et al., 2008; Torres-Romero and Olalla-Tárraga, 2015), whereas others have highlighted a negative relationship between species richness and remaining habitat (Da Silva et al., 2018; Velazco et al., 2023, 2019). The presence of large regions with high species richness and high remaining habitat in Colombia suggests that many areas still hold potential for conservation efforts to protect snakes and other organism groups.

Priority areas for conservation are concentrated in the Pacific, Amazon, Orinoco, and northern Caribbean regions. Although the overlap with PA was low, it increased with the inclusion of the OECMs (Fig. S5). The largest number of communities and indigenous lands in the country are found in the northern Caribbean, Orinoco, and Amazon regions (DANE, 2012), in addition, the Pacific, Amazon, and Orinoco are regions with the highest proportion of natural cover (Fundación Gaia Amazonas, 2023). These regions represent an opportunity to identify new OECMs that, together with governmental actions will contribute to consolidating more ambitious conservation goals (Shiono et al., 2021). Biodiversity and cultural diversity are interconnected, and healthy ecosystems are the basis for the existence of indigenous peoples and local communities (Levis et al., 2024). Therefore, identifying new OECMs may represent an opportunity to transform conservation policies and practices and to recognize the contributions of indigenous peoples and local communities. (Jonas et al., 2017; Levis et al., 2024). OECMs are changing the paradigm regarding conserving biodiversity, as they generate more inclusive and representative systems that evidence multiple strategies, actors, and institutional governance arrangements at the local scale (Jonas et al., 2014).