Colombia is located in the northwest part of South America with an area of 1.1 million km2. The country is divided into five main biogeographic regions: The Andean, Caribbean, Pacific, Colombian Amazon, and Orinoquia plains (Etter et al., 2006).

Due to its complex geologic history regarding the uplifting of the Andes Mountain range, and its tropical weather, Colombia has a high ecosystem diversity, including almost all biomes described for the planet (Etter, 1993). To protect and conserve its biodiversity, Colombia has created a National System of Protected Areas (SINAP) with around 1000 protected areas grouped into 14 protection categories; this represents around 16% of the terrestrial land of the territory (RUNAP, n.d.). Colombian forests cover 52% of the territory (IDEAM et al., 2021) and, despite their degree of protection, continue to be lost and degraded (even within PAs) due to their conversion to pastures for cattle ranching, the expansion of the agricultural frontier, as well as the cultivation of illegal crops (Armenteras et al., 2013; Etter et al., 2006).

The selection of species for large-scale connectivity analysis is a complex task. Given the national scale of our study, we adopted a multispecies perspective and therefore selected species that jointly represented or encompassed the ecological needs of many others at country level (sensu Lambeck, 1997). We further grouped species together to create “species ecological profiles” (de la Fuente et al., 2018; Gurrutxaga et al., 2010; Vos et al., 2001) based on three criteria: (1) preferred land cover, which constitutes the species residential habitat, (2) altitudinal distribution ranges of the species, and (3) dispersal distances.

Connectivity within and among PAs is threatened by forest changes and habitat loss. Therefore, our analysis included species with restricted distributions as well as those highly sensitive to habitat fragmentation (Henle et al., 2004). The initial step in species selection was the compilation of the list of threatened species of Colombia (Colombian Resolution 1912 of 2017). Mammals were selected from this list as our focal taxon. They are often used in connectivity studies because of their susceptibility to linear infrastructure and the potential for other species to use the corridors designed for them (Gurrutxaga et al., 2011).

From this initial list, we excluded species with wide distributions (i.e., generalists) from the analysis because our focus was on species with specific habitat selection (forest-dependent species). Additionally, we excluded species with small body mass because their dispersal capabilities are smaller than the grain and scale considered in this study (1 km2). The data on habitat selection and dispersal capabilities of the species were obtained from relevant scientific literature (see Appendix A Table A1). After applying these filters, we selected 16 threatened mammal species and grouped them into 4 ecological profiles (see Table 1). Individual home range values were averaged to obtain the mean home range by ecoprofile. Subsequently, this value was used to calculate the maximum dispersal distance (see Table 1).

To assign the species to a specific ecological profile, we considered two factors: (1) the maximum altitude value they occupy as for their specific altitudinal range, and (2) the preferred land cover, following the classification by Rodríguez-Eraso et al. (2006) where Andean ecosystems are reviewed and described. The term “sub-Andean forest” was used by Cuatrecasas (1958) to define the humid forests between 1000 and 2400 m. Furthermore, Cuatrecasas (1939) described the limits of the Andean forests as being between 2400 and 3000 m, and the high Andean between 3000 and 3300 m. The resulting ecological profiles are presented in Table 1.

Habitat patches are defined as adjacent areas that meet specific habitat requirements for a target species; they represent unique zones where species can reproduce, find prey supply, and shelter (Wade et al., 2015). We identified potential habitat patches for our connectivity analyses by selecting the patches of preferred land cover and the targeted altitudinal ranges for each ecological profile. The patches of preferred land cover were obtained from the national map of land cover types for Colombia (IDEAM, 2021) at a scale of 1:100,000. To facilitate the overall analysis, we reclassified the original 54 categories of land cover into 13 different classes. Elevation data was derived from satellite images obtained from the NASA Shuttle Radar Topography Mission (NASA JPL, 2013) with a resolution of 30 m. We reclassified the elevation data into 7 categories. Altitudinal ranges for each species were obtained from the database “Mamíferos de Colombia, version 1.12” (Ramírez-Chaves et al., 2021). Subsequently, we overlapped the resulting potential habitat patches with the National System of Protected Areas for Colombia (SINAP).

Finally, we defined nodes inside these habitat patches within PAs as the starting and ending points of potential movements of individuals of the target species. As a general rule, we placed one node in the centroid of each habitat patch. However, in the study area there are some protected areas and patches disproportionally big in relation to the others and placing just one node in their centroid could bias the spatial definition of the movement routes. Therefore, we created a grid with a cell size of 50 × 50 km and intersected it with the selected habitat patches. In these specific cases, we considered intrapatch connectivity in our analyses (sensu Saura and Rubio, 2010) by assuming each grid cell within the large forest patch to represent one independent node for subsequent analyses.

For the generation of resistance surfaces, we considered two environmental variables: land cover type and elevation. These layers were parameterized by assigning a cost or resistance value to each 1 × 1 km pixel across the study area, where low resistance represents ease of movement, and high values denote restriction to movement. Empirical data for the species, such as presence-absence points, capture pathways, or genetic data are lacking, therefore, a one-stage expert’s consultation approach was used, where the expert’s opinion is considered to derive the final resistance surface (Zeller et al., 2012).

A group of 22 experts on Colombian mammals was contacted to participate in the parameterization of resistance surfaces. To collect the expert’s opinion, we created a questionnaire with the categories of land cover type and altitudinal ranges. We asked species’ experts to assign resistance values to each combination of species-landcover and species-altitude categories. Resistance values ranged from 1 to 100, with 1 representing the habitat for the species and 100 the maximum resistance to movement; it was also possible to assign a value of 1000 to those categories that the experts considered insurmountable barriers to movement (e.g., the Andean range for lowland species). The results obtained from the questionnaires were averaged for land cover type and altitude categories within each ecological profile.

For this study, we assumed that species select residence habitat with the same criteria as dispersal habitat. Furthermore, it is important to notice that although resistance values were assessed individually for each species, these values are intended to represent resistance to movement for a group of species, in this case an ecological profile, and do not pretend to model the dispersal behavior of a single species.

To model potential movement routes among habitat patches, we first applied a Least Cost Path (LCP) analysis using the Linkage Mapper toolbox (McRae and Kavanagh, 2011 ), and then assessed the priority of each resulting corridor within a conservation scenario by using the software Conefor 2.6 (Saura and Torné, 2009).

In the LCP analysis, adjacent nodes were identified based on both cost of movement (given only by the resistance layer of each ecoprofile) and shortest euclidean distance. If two nodes were not adjacent by cost-weighted distance, the shortest euclidean distance was considered. It is important to clarify that, for this step, we did not set a limit to the cost of movement or distance when we ran the analysis, because we wanted to construct a thorough network identifying all the potential corridors connecting adjacent habitat patches.

We refer to the LCPs as corridors, although we acknowledge that these lines would be more accurately regarded as the central axes of the swaths of land around which movements are more likely to occur. After identifying all the corridors between adjacent nodes, we assessed the priority of each within a conservation scenario. This assessment was based on their significance for maintaining habitat availability at the landscape scale, considering the study area as a whole. For this prioritization we applied the Probability of Connectivity index (PC) (Pascual-Hortal and Saura, 2006; Saura and Pascual-Hortal, 2007; Saura and Rubio, 2010). Each corridor was ranked within the corresponding ecoprofile according to the decrease of the PC index value after its removal (dPC) (see Box 1).

To parameterize the dispersal capability in the Conefor calculations, for each ecoprofile we took the maximum dispersal distance across all its species. This distance can be regarded as a dispersal event of exceedingly low probability and accordingly we set its probability to 0.0001 (see Box 1, equation 2). We estimated the dispersal distance of each ecoprofile using the equation proposed for mammals by Bowman et al. (2002). Subsequently, we transformed these distances into cost units by multiplying the mean cost value extracted from the corresponding resistance surface (see Cost Weighted Distance in Appendix A Table A3).

We selected as priority corridors for each EP the 30% of the LCPs with the highest dPC values, assuming these could potentially constitute the backbone for a national connectivity strategy. Then, we assessed the landscape context around these corridors. For that, we created a buffer of 2 km around all of them (following Beier, 2019), namely corridor swaths, and extracted data about the land cover types and potential human disturbances through the Human Footprint Index (HFI) (Correa Ayram et al., 2020) within these swaths.

The habitat patches identified for the lowland and sub-Andean ecological profiles (EP1 and EP2) were distributed across all five Colombian biogeographical regions, with larger habitat patches located in the Amazon (see Fig. 1A and B). The smallest mean habitat patch size was identified for Andean and high Andean small-medium dispersal ecoprofile EP3 (12,000 ha), for which we only identified a few isolated suitable habitat patches in regions other than the Andes (Fig. 1C). In the Andean region the EP4 habitat patches were mainly distributed along the eastern and western slopes of the mountain range, and the Caribbean nodes were inside the NNP Sierra Nevada de Santa Marta (Fig. 1D). The figures and summary tables for each ecoprofile can be found in Appendix A Figure A1 and Table A2.

Fig. 1.

Least Cost Paths (LCPs), calculated based on resistance layers for each ecoprofile, linking adjacent nodes (gray dots) of forest habitat patches (green polygons). Gray lines represent the LCPs that surpassed the maximum dispersal distance of the ecoprofile (see Appendix A Table A3) while yellow lines indicate those potential corridors where Cost Weighted Distance (CWD) and Euclidean distance do not exceed those maximum dispersal values. A. Ecoprofile 1: Lowland and sub-Andean small and medium dispersal capability; B. Ecoprofile 2: Lowland and sub-Andean large dispersal capability; C. Ecoprofile 3: Andean and high Andean small and medium dispersal capability; D. Ecoprofile 4: Andean and high Andean large dispersal capability.

(1.32MB).

The experts’ questionnaires had a return rate of 32% (7 experts out of 22). In general, for all ecoprofiles, experts predicted high resistance for urban areas, zones without vegetation, agiculture, agricultural mosaics and pastures (see values in Appendix B). Furthermore, the land cover plantations had high values especially for EP1 and EP3. Water land covers were all set with a value of 1000 as insurmountable barriers.

Other non-forest land covers, like grasslands and shrublands as well as secondary vegetation, were parameterized with values around 30 and 50, which allows for medium permeability in areas with these types of land covers. Evidently, EPs 1 and 2 had lower resistance for lowland areas and the opposite was true for EPs 3 and 4, as shown in Appendix A Figure A2.

Our results showed an extensive network of potential corridors between selected nodes for each ecoprofile, which had average lengths ranging from 70 to 100 km (see Fig. 1). A great number of potential corridors were identified inside large habitat patches, especially for lowland and sub-Andean ecoprofiles (1586 and 942 corridors, respectively; Fig. 1A and B). The ecoprofile EP2 showed the highest average length of corridors among all, with values around 100 km. Most corridors for EP3 (1332 in total) concentrated over the Andean region, however, our analysis also showed a few isolated potential corridors linking habitat patches from the Orinoquia region to the Andes (see Fig. 1C); these potential corridors are about 500 km long and it would be virtually impossible for the species of this ecoprofile to use them, since they have small and medium dispersal capabilities and there are not habitat patches that could act as stepping stones between the nodes that these corridors link. Lastly, the ecoprofile EP4 presented the lowest number of potential corridors (571 in total) with an average length of 79 km, these are located mainly in the western and eastern slopes of the Andean Mountain range (Fig. 1D).

We identified a large proportion of priority corridors for EP1 (as ranked by dPC values) within the Amazon, Orinoquia and Pacific regions (53%, 55%, 83% of the region’s total priority corridors, respectively). Some of these corridors were concentrated inside large habitat patches, especially those located within PAs of the Amazon region, such as: National Natural Park (NNP) La Paya (A), NNP Serranía de Chiribiquete (B), NNP Sierra de la Macarena, Tinigua, Picachos (C), and Natural reserve Nukak (D) (see Fig. 2). Similar results were obtained for species grouped in EP2. The priority corridors of ecoprofile EP3 were mainly located in the Andean, Caribbean and Pacific regions, showing the lowest average length among all ecoprofiles (around 43 km). Lastly, we identified several priority corridors for ecoprofile EP4 connecting habitat patches within the Andean and Caribbean region. For the latter, the corridors were located especially inside NNP Sierra Nevada de Santa Marta and the eastern slope of the Andean region (see Northeast area Fig. 2).

Fig. 2.

Priority corridors for conservation represented over the values of the Human Footprint Index (Correa Ayram et al., 2020). Colored lines represent the central axis of the potential priority corridors obtained from the prioritization analysis in Conefor (top 30% dPC values for each ecoprofile; see details in Box 1). The colored doughnuts show the percentage of the 2-km corridor’s swaths within each biogeographic region that corresponds to each ecological profile. Labels: A. National Natural Park (NNP) La Paya; B. NNP Serranía de Chiribiquete; C. NNP Picachos, NNP Tinigua, NNP Sierra de la Macarena; D. Natural reserve Nukak.

(0.69MB).

Overall, all priority corridors were located in areas with low HFI, showing mean values ranging from 1 to 20. The highest values of HFI were concentrated in the Andean region (see Fig. 3A), especially for Andean species with small, medium, and large dispersal capabilities (EP3 and EP4).

Fig. 3.

(A) Human Footprint Index (HFI) (Correa Ayram et al., 2020) values within the 2-km priority corridor's swaths identified to connect protected areas, detailed by ecoprofile (EP) and biogeographic region (mean values as circles and standard deviation as lines). (B) Percentage of Land Cover type (IDEAM, 2021) within 2-km corridor’s swaths by ecoprofile (EP).

(0.44MB).

The corridors within the Caribbean region also exhibited HFI values similar to the corridors located in the Andes, except for EP2 which showed a much lower value of HFI (see Fig. 2A). On the other hand, the corridors located in the Amazon and Orinoquia region traversed through areas with very low values of human pressure (mean HFI of 2–3).

The priority corridors swaths traversed mainly through dense forests, with percentages of area above 90 % for lowland and sub-Andean species (EP1 and EP2), and above 50 % for Andean and high Andean species (EP3 and EP4). For the latter, priority corridors traversed through other land cover types such as agriculture mosaics, grasslands, secondary vegetation, and pastures, with values around 20% of the total area (see Fig. 3B).