Omnidirectional connectivity for the Andean bear (Tremarctos ornatus) across the Colombian Andes

Conserving or restoring connectivity is a common objective of landscape-scale conservation initiatives. However, precise species occurrence or movement data to inform or validate spatial models are often lacking. Our objectives were to (1) produce the first approximation of country-wide connectivity for Andean bears (Tremarctos ornatus) in Colombia and (2) demonstrate a novel approach for model validation which uses publicly available web and social media records of a flagship species. We used general knowledge about Andean bear habitat associations and indices of ecological integrity to construct a resistance surface across the Colombian Andes. We used this resistance surface to model omnidirectional connectivity using circuit theory. We validated our model with coarse location data acquired from local news stories and social media posts. Our model was most sensitive to changes in the resistance values of agricultural landcover and the mid-elevational zone, but uncertainty analysis demonstrated these changes had little impact on our conclusions regarding the municipalities most conducive to Andean bear movement. Just over one-third of those areas most conducive to Andean bear movement were within protected areas, while 8% coincided with agricultural landcover. We constructed a model of connectivity that did not rely on independent, empirically derived location data. Our model is coarse (1 km resolution) but can still provide useful information to practitioners in Colombia who are working with scarce ecological data. More information about how Andean bears move through agricultural landscapes would help improve our understanding of connectivity for this species in Colombia.


Introduction
Conservation practitioners have frequently identified the improvement or maintenance of landscape connectivity as an overarching objective in landscapescale conservation initiatives (Koen et al. 2019;e.g., IUCN 2020). These initiatives rely on a variety of spatial modeling approaches to determine where conservation efforts would be most effective or advantageous. The most popular models for prioritizing landscape elements for connectivity conservation are premised on cost-based analyses (Cushman et al. 2013). These approaches require the analyst to identify specific landscape characteristics (e.g., landcover, elevation, road density) that enable or hinder movement to differing degrees. These relative differences in movement conduciveness are parameterized into connectivity models as resistance values or ''costs'' to movement; landcovers considered permeable to movement have low resistance, while impermeable ones have high resistance. Most early connectivity models were least-cost path analyses (LCPs) which identified linear elements in the landscape which were most conducive to movement and provided the literal ''path of least resistance'' between two focal habitat patches.
More recently, ecologists have turned to connectivity models guided by electrical circuit theory (McRae et al. 2008). These models, usually operationalized within the open-source program, Circuitscape Dickson et al. 2019), take advantage of parallels between electrical current and animal movement (or gene flow) across a landscape (McRae et al. 2008). Contrary to LCPs, circuitbased models acknowledge and account for multiple dispersal pathways that contribute to functional connectivity across a given landscape in the same way that multiple wires in an electrical network simultaneously contribute to lessen the effective resistance between nodes (McRae et al. 2008). While LCPs assume that individual animals have perfect knowledge of the landscape, circuit-based models draw on correlated random walk theory and thus assume no prior knowledge . Instead, each movement is determined by what the theoretical organism is confronted with in its immediate vicinity, though a general directionality is maintained. Areas identified by circuit-based models as having the highest current are those which have the highest probability of use by ''random walkers'' (McRae et al. 2008). Dickson et al. (2019) reviewed publications citing Circuitscape and found that the method has been used for connectivity analyses on every continent. Circuit-based models have successfully predicted locations of Canada lynx (Lynx canadensis; Walpole et al. 2012), ovenbirds (Seiurus aurocapilla;St-Louis et al. 2014), African wild dogs (Lycaon pictus), and cheetahs (Acinonyx jubatus; Jackson et al. 2016) and have been validated by numerous genetic studies (e.g., Blair and Melnick 2012;Devitt et al. 2013). Circuit-based models have outperformed LCPs in predicting functional connectivity for moose (Alces alces; Laliberté and St-Laurent 2020), dispersal by wolverines (Gulo gulo;McClure et al. 2016), and road-crossing by black bears (Ursus americanus; Zeller et al. 2020).
Circuit-based models generally begin with the delineation of source and target patches or ''nodes'' (e.g., two national parks between which conservation practitioners would like to improve functional connectivity), but in 2014, Koen et al. introduced an approach that does not require focal source and target patches. Instead their approach incorporates many randomized nodes located beyond the perimeter of the study area. This approach creates a model of omnidirectional connectivity within a landscape (some researchers refer to this as ''wall-to-wall'' connectivity) (see also Pelletier et al. 2014 for an alternate approach to modeling omnidirectional connectivity). Omnidirectional modeling is ideal for scenarios in which there are no a priori justifications for source and target patches (e.g., if a practitioner wishes to assess connectivity across the entire state of California) and minimizes the need for potentially complicated siteselection processes (see Beier et al. 2011) that have a disproportionately high impact on the final model. Finally, the approach of placing nodes in a buffer beyond the perimeter of the study area resolves issues of node-placement bias (i.e., the propensity for current values to be artificially inflated near source and target nodes) (Koen et al. 2014).
Practitioners working towards the conservation or restoration of landscape connectivity are often challenged by the need for empirical data required to both develop and then validate models of connectivity. For example, setting resistance values is considered to be ''probably the most important bottle-neck for applying [cost-based analyses]'' (Adriaenson et al. 2003). Movement data (sometimes referred to as pathway data) (e.g., from GPS collars) is generally considered to be the gold-standard for focal-species based connectivity models. However, studies that seek movement data are notoriously expensive and even highquality occurrence data are limited for most species. Andean bears (Tremarctos ornatus)-the only extant ursid in South America-are one such example of a data-poor species Falconi et al. 2020). This is particularly true within the country of Colombia where many regional environmental authorities have neither the time nor the financial resources to carry out research on Andean bears yet are responsible for implementing conservation programs for the species within their jurisdictions (Velásquez Durán 2018). However, Koen et al. (2014) demonstrated that even a generic resistance surface for species that preferred ''natural'' landcover types resulted in a Circuitscape model adept at predicting fisher (Pekania pennanti) habitat use and amphibian and reptile road mortality in eastern Ontario. This type of ''naturalness'' model is sometimes referred to as a connectivity model of ecological integrity (e.g., Beier et al. 2011). Connectivity models of ecological integrity have been shown to be particularly suited to modeling landscape connectivity for large-bodied, vagile species (Krosby et al. 2015). Thus, we saw the Koen et al. (2014) approach to modeling omnidirectional connectivity with ecological integrity as a potentially useful approach for creating a countrywide connectivity model for Andean bears in Colombia despite data limitations.
Model validation remains the final obstacle to creating reliable models of connectivity for focal species and ecological integrity alike. Indeed, Laliberte and St-Laurent (2020) very recently wrote of the proliferation of connectivity models that were not empirically validated, calling model validation the ''Achilles' heel of landscape connectivity mapping.'' However, we argue that broad-scale connectivity models do not necessarily need fine-scale data for validation. The model only needs to be accurate to the level at which management actions or decisions would be applied. Thus, we searched for unconventional forms of data that could be used to validate our model at the level of the municipality in Colombia. These data were found in the form of publicly available records on the web, from news media, and from various social media platforms. Andean bears are classified by the International Union for the Conservation of Nature as a species vulnerable to extinction (Velez-Liendo and Garcia-Rangel 2017). An estimated 3000-6000 Andean bears occur in Colombia (Ruiz-García 2003). In 2001, the Ministry of the Environment in Colombia commissioned a panel of experts to assemble the National Program for the Conservation of the Andean bear (Mayr Maldonado 2001). Various state and private entities alike have worked towards the objectives outlined within this program over the last two decades, though implementation of the national program has been slow (Rodríguez-Castro et al. 2015). Many practitioners have spoken of the challenges of implementing conservation efforts for a species for which they have such scarce ecological data (pers. observ.). A model of landscape connectivity helps to address this issue of data scarcity by, for example, elucidating those areas in the landscape that if protected or restored could enhance the population viability of the species across this fragmented landscape (Kattan et al. 2004).
We created a coarse (1 km resolution), broad-scale connectivity model for Andean bears based on general knowledge about their preferences for high-elevation habitats and known avoidance of human-dominated landscapes. We had two main objectives guiding this research: (1) produce the first approximation of country-wide connectivity for Andean bears in Colombia and (2) demonstrate a novel approach for model validation which relies on publicly available web records and social media data of a flagship species. In addition to these two objectives, we also sought to identify any challenges or disparate outcomes when applying the Koen et al. (2014) approach to a much broader spatial extent (our study area encompassed over 400,000 square km). Though their methodology had been cited over 100 times at time of writing, very few researchers conducted their own independent analyses of requisite buffer widths and number of nodes in addressing node-placement bias. Thus, we do not yet know if their recommendations are consistent for mapping omnidirectional connectivity in study areas covering larger spatial extents. Though our model was developed for Andean bears, it also serves as an approximation of the connectivity of ecological integrity across the Colombian Andes. This model could serve as a valuable decision-support tool for practitioners working with scarce ecological data in the country. Moreover, the approach demonstrated herein is extremely cost-effective and can serve as a first step towards identifying priority conservation areas in other data-scarce regions and with other charismatic species.

Colombian Andes
The Colombian Andes have three distinct mountain ranges (or cordilleras) that run from the southwest to the northeastern borders of the country. The Colombian Andes collectively account for approximately one-fourth of the country's total surface area and are acknowledged to be a biodiversity hotspot and thus a high priority for conservation action (Mittermeier et al. 1999;Armenteras et al. 2003). The Colombian Andes contain a diverse array of tropical ecosystems including cloud forest, wetlands, and high-altitude shrub ecosystems known as pa´ramo.
Andean ecosystems have experienced extensive landcover change since Colombia's colonization by Spain in the early 1500s, the majority of which was driven by the introduction and rapid expansion of cattle ranching (Etter et al. 2008). Exponential human population growth in the latter half of the twentieth century and rural-to-urban migration are also important factors in recent landcover change (Etter et al. 2008). Over 60% of Colombia's 50 million residents currently reside in the Andes and the intermountain valleys. Andean cloud forests are estimated to cover less than 50% of their historic extent while 15% of Colombian pa´ramo has been converted to agricultural uses (Llambí et al. 2019). The cordilleras vary in their degree of landcover conversion, with the eastern cordillera being the most developed of the three. Climate change further threatens these high-altitude ecosystems (Urrutia and Vuille 2009;Buytaert et al. 2011).

The Andean bear
Andean bears, also known as spectacled bears (in Spanish, osos andinos or osos de anteojos), are the last surviving lineage of the subfamily Tremarctinae (García-Rangel 2012). They are found across the northern Andes in Venezuela, Colombia, Peru, Ecuador, Bolivia, and possibly Argentina (Cosse et al. 2014;Velez-Liendo and García-Rangel 2017). In Colombia the species primarily occupies Andean cloud forest and pa´ramo (Peyton 1999). However, research conducted in other countries documented this species in several other ecosystems such as tropical dry forest (Peru; Kleiner et al. 2018), steppe grasslands (Peru; Peyton 1987), and elfin forest (Bolivia; Ríos-Uzeda et al. 2006). The lowest elevation at which Andean bears have been recorded was 200 m (Velez-Liendo and García-Rangel 2017), but in Colombia as elsewhere they are most often found in mountainous ecosystems above 1200 m (Peyton 1999;Ríos-Uzeda 2006).
Andean bears have experienced extensive habitat loss and fragmentation across their range. Kattan et al. (2004) estimated that remaining habitat for Andean bears comprises approximately 200,000 km 2 -42% of the suspected original extent of the species. Kattan et al. (2004) further found that this habitat was fragmented into 113 distinct patches, with most (56%) intact patches smaller than 500 km 2 (''intact'' indicating they were not intersected by roads). The minimum patch size believed to support viable Andean bear populations is between 1200 and 1900 km 2 (Yerena 1998;Peyton 1999). In addition to issues of habitat loss, Andean bears are further threatened by mortality due to human-wildlife conflict (Velez-Liendo and Garcia-Rangel 2017). Though Andean bears are largely herbivorous, they do scavenge, and some individuals will occasionally kill cattle (Parra-Romero et al. 2019). They can also cause significant damage to crops (Peyton 1980;Escobar-Lasso et al. 2020). Thus, many campesinos (people living in rural areas) across the Andes see Andean bears as a threat to their livelihoods and persecute them as such (Jorgenson and Sandoval-A 2005;Goldstein et al. 2006;Parra-Romero 2011;Figueroa 2015).

Delineating the study area
We acquired broad-scale data on the range of Andean bears in Colombia from the International Union for the Conservation of Nature (IUCN 2017). This dataset consisted of polygons delineating known or suspected Andean bear presence across all three cordilleras of the Colombian Andes. We cross-validated these data with interviews conducted with conservation practitioners in Colombia during 2018-2019 (Hohbein 2021). We subsequently added to this dataset polygons representing Serrania del Perija Regional Natural Park and Paramillo National Natural Park which were known by local practitioners to have Andean bears, but which were absent from the IUCN dataset (presence in Serrania del Perija was also described in Rodríguez et al. 2019). We used the ArcMap tool ''Minimum Bounding Geometry'' (convex hull) to incorporate all range polygons into a single feature. The resulting polygon was then clipped to the extent of Colombia to function as our study area. Our final study area encompassed 417,895 km 2 .
Calculating the resistance surface

Variables included
We used landcover and transportation corridors as our two indicators of anthropogenic impact. We also included elevation as a consideration for connectivity as Andean bears are known to primarily occur in highelevation ecosystems. We combined rasters representing the resistances from each of these variables to create a cumulative resistance surface layer for Andean bears within Colombia. More information on our process for assigning resistance values and creating each of the resistance rasters is provided below.

Choosing resistance values
True resistances to movement for organisms are rarely (if ever) known. Methods for interpreting suspected resistances to movement vary substantially and require numerous-and often subjective-decisions by the analyst (Beier et al. 2008). For example, should the costs to movement range from 0 to 1 or from 1 to 1000? Is agricultural land twice as hard to move through than natural landcover or five times harder? Fortunately, Bowman et al. (2020) found that so long as the relative ranks of various landscape variables are accurate, circuit-based models are relatively resilient to changes in cost assignments. We chose to follow Koen et al. (2014) in their assignments of resistance values, with the lowest possible resistance assignment set to 10, moderate resistance set to 100, and high resistance to movement set to 1000. However, we felt we needed one more degree of flexibility for resistance assignments and added an option for a resistance value of 500 to represent moderately high resistance. We did not include any barriers to movement in our analysis.

Landcover
We derived our landcover resistance surface from the most recent classified landcover dataset available from Colombia's Institute of Hydrology, Meteorology, and Environmental Studies (Instituto de Hidrologı´a Mete-orologı´a y Estudios Ambientales [IDEAM] 2014). This dataset consisted of 53 landcover classes at a resolution of 1:100,000 from 2010 to 2012. We grouped these original classes into 6 broad categories and assigned these as having one of the four possible resistance values: natural (10), agriculture (100), inland water (including rivers [ 50 m across and lakes; 500), rural (500), urban areas (1000), and other ''unnatural'' (1000) (see Online Resource 1 for full list of classifications). We then converted this feature class to a raster layer with 1000 m resolution; cells were assigned to the landcover class that had the largest area within the cell. We created a distinct river raster layer prior to conversion to better maintain feature continuity and address issues of diagonal discontinuities known to be problematic for eight-neighbor algorithms (Adriaensen et al. 2003). This river raster (also given a 1000 m resolution) was then integrated with the other comprehensive landcover raster layer with the river raster given priority.

Elevation
Our elevation data was comprised of a 30 arc-second digital elevation model of South America from the U.S. Geological Survey's Center for Earth Resources and Observation and Science (EROS 1996). Because the lowest known Andean bear sighting was at 200 m, we classified anything below this value as having high resistance (1000). Elevations between 200 and 1200 m were assigned medium resistance (100), while elevations above 1200 m (where Andean bears are most commonly observed) were assigned a low resistance value (10). Though Andean bears are not often observed at elevations above 4000 m, we decided not to assign higher resistance values to these higher elevations because these higher elevation areas are likely to be important for future connectivity given predicted climate change impacts on mountainous ecosystems (Walther et al. 2002;Urrutia and Vuille 2009).

Transportation corridors
We acquired primary Colombian transportation routes from the Colombian National Institute of Roads (Instituto Nacional de Vias [INVIAS] 2019) which we converted to a raster layer with a 1000 m resolution. As with rivers, we addressed issues of diagonal discontinuities to ensure there were no ''holes'' in these barriers (Adriaensen et al. 2003). This dataset included only those major roads that connected larger Colombian cities (i.e., primary roads only); thus, we assumed at least a moderate volume of traffic. We assigned all roads included in this data set as having a resistance of 500.

Resistance values beyond Colombia
Arbitrary jurisdictional and study area boundaries can create artificial effects in Circuitscape models (Koen et al. 2010(Koen et al. , 2014. We needed both an ocean resistance layer beyond Colombia's shoreline and a suitable substitute for resistances beyond the terrestrial borders of Colombia for which we did not have data. To create the ocean resistance layer, we assigned oceans within a 100-km buffer around Colombia a resistance value of 2501 (to exceed the highest cumulative resistance of 2500) (see Table 1 for full list of variables and corresponding resistance values). Koen et al. (2010) determined that random resistance layers placed beyond the limits of the study area were as effective as resistance surfaces derived from real data in eliminating artificial boundary effects on resultant connectivity models. Thus, we followed their methodology in creating a random resistance layer that would substitute for unknown resistance values in the neighboring countries of Panama, Venezuela, and Ecuador. To do this, we used the random.raster function from the spatialEco package (v. 1.3-1, Evans 2020) in R (v. 1.1.383, R Core Team 2019) to create a raster layer that extended 100 km beyond Colombia's terrestrial borders. Resistance values for this new random raster were sampled from the cumulative resistance surface within our delineated study area. The resultant random raster and ocean raster were joined to the cumulative resistance layer for Colombia to create our final resistance surface (Fig. 1).

Comparing Circuitscape software and Julia update
In 2019, the Circuitscape connectivity program was upgraded to operate with the Julia programming language, enhancing processing power and allowing for applications to much larger spatial extents and/or with higher resolution data ). To ensure model agreement as well as eliminate the possibility for user error in Julia script input, we performed a test on a portion of our study area (84,912 km 2 ) before initiating the buffer width and node analysis. For this test, we calculated the cumulative current map for the smaller area using the original Circuitscape software and compared this against the cumulative current map produced by the Julia upgrade  to ensure an exact match in current density estimates.
Determining appropriate buffer width and number of nodes Following Koen et al. (2014), we sought to remove the effects of node-placement bias by randomly assigning node locations beyond the extent of the study area and eliminating the buffer which contained these effects before assessing landscape connectivity. We used a 3-step approach to determine the appropriate buffer size and number of nodes: (1) determine a suitably large buffer width using an arbitrary number of nodes (n = 50); (2) determine the optimal number of nodes at the identified buffer width; and then (3) repeat the buffer analysis with the corrected number of nodes.

Buffer width
Following the recommendation of Koen et al. (2014), we chose to start our buffer analysis with 50 nodes (Koen et al. considered 50 nodes to be their ''full pairwise map''). We tested buffer widths between 10 and 100 km (approximately 2% and 20% of maximum study area width, respectively) at 10 km intervals. For each buffer width tested, we randomly generated 50 nodes along the outer perimeter of the buffer; we maintained a minimum of 10 km between neighboring nodes. We used Circuitscape in Julia to connect all possible node pairs using the eight-neighbor rule, imported the resultant cumulative current map into ArcGIS, and removed the buffers. This allowed us to compare landscape connectivity only within our delineated study area and with a reduced bias from node placement following buffer removal (Koen et al. 2014). The cumulative current map produced using each buffer width was compared against the cumulative current map produced using the largest buffer size tested (100 km) with Pearson's correlation coefficient (r). To determine whether the 100-km buffer was sufficiently wide to have removed biases from node placement, we looked for evidence of negligible or no improvement in Pearson's r following an increase in buffer width. To do this, we first calculated the rate of change (m) in Pearson's r between all consecutive pairs of current maps produced at all buffer widths tested. For example, the rate of change between  We then calculated one-sided 3-point moving window averages (MWA) for the slope values. For example, to calculate the one-sided 3-point MWA for the 70-km buffer: When the 3-point MWA fell below 0.001 for two values in a row, we considered this to be sufficient evidence that further increasing the buffer width would result in negligible improvement to Pearson's r. Had this not occurred during our tests, we would have needed to repeat the process with an even larger maximum buffer (e.g., 200 km). However, this procedure demonstrated that the 100-km buffer was sufficiently wide (the 3-point MWA fell below 0.001 a second time at 70 km).

Number of nodes
Similar to the procedure above, we tested the effects of increasing numbers of nodes on Pearson's r. To do this, we first created 100 random nodes around the outer perimeter of a 100-km buffer. We then randomly selected increasingly larger subsets of nodes from this pool of candidates for each consecutive model. We initially began our testing with 30 nodes, and we increased the number of nodes by 10 until we reached the full set of 100 nodes. We compared all cumulative current maps produced with differing numbers of nodes against the cumulative current map produced with the full set of 100 nodes. We detected a leveling off in the Pearson's correlation coefficients between 70 and 80 nodes. Following this result, we then tested between 60 and 90 nodes by increments of 2. We followed a similar procedure as used above in the buffer analysis (one-sided 3-point MWA of slope values) for identifying the optimal number of nodes (n = 78; 3003 unique pairs) (Fig. 2a).

Confirming buffer width selection
Following the identification of a suitable number of nodes for this analysis, we reran the buffer analysis to determine whether the appropriate buffer size would change given the revised number of nodes used. We reran Circuitscape in Julia with 78 nodes at 5-km buffer width intervals that ranged between 10 and 100 km; for each buffer width tested, we generated 78 new nodes along the perimeter of the buffer that were at least 10 km apart from neighbors. We again calculated Pearson's r for the current maps produced for the study area using all buffer widths compared against the current map produced with the largest buffer tested: 100 km. We calculated one-sided 3-point MWAs for all consecutive rates of change in Pearson's r. We determined the buffer was sufficiently wide to have removed node-placement bias at 45-km, the point at which two consecutive MWAs fell below 0.001 (Fig. 2b). Thus, the final current map was produced using 78 nodes and a 45-km buffer.

Validating model
Andean bears are theoretically more likely to traverse through those areas in our model with higher current densities than they would those areas with lower current densities. Therefore, we would expect Andean bear presence and movement data to be associated with locales with higher mean current. However, data regarding precise Andean bear movements or even presences and absences are extremely limited. Not having access to such data throughout the study area, we validated our connectivity model with a novel approach that allowed us to gain coarse location information using publicly available web records. We searched the web for various mentions of Andean bear sightings, poaching reports, or publicly accessible location data. In November 2019, we searched GoogleÒ, TwitterÒ, InstagramÒ, and FacebookÒ for the terms ' We examined all web results, video results, and tweets for mentions of specific municipalities or neighborhoods where Andean bears had either been seen or recently poached. Additionally, we searched Colombia's Biodiversity Information System (Sistema de Informacio´n sobre Biodiversidad de Colombia [SiB]) and the Global Biodiversity Information Facility (GBIF) for publicly available location data for Andean bears. We included all referenced locations irrespective of the publication or observation date. Through these combined processes, we identified 421 references to places with Andean bears within 110 different records (numerous place references often occurred within a single document, webpage, tweet, etc.). Of these 421 references, 235 were coordinates from across SiB and GBIF (GBIF 2021). The remainder of these records were from numerous different sources, including agency reports and websites, local news outlets, and social media (see Online Resource 2 for citations). We further subsidized these 421 references with an additional 106 references made by conservation practitioners during interviews conducted in 2018 and 2019 (Hohbein 2021). Cumulatively, these 527 references represented 155 different municipalities in Colombia.
We calculated the mean current value from our final connectivity map across the 155 municipalities with records of Andean bear presence and the mean current of all other municipalities within our delineated study area for which we found no records (n = 775). An F-test on the variances of each group of municipalities showed support for equal variance within each group. Following confirmation of equal variance, we performed a one-sided T-test to assess whether mean current was higher within those municipalities with Andean bear sightings than those without such records. We repeated this process for those SiB and GBIF records which included geographic coordinates by calculating mean current within a 1-km buffer of each record (n = 235) and comparing this against the mean current within a 1-km buffer of randomly generated points (n = 1175) distributed throughout the study area.

Sensitivity analysis
We conducted a sensitivity analysis to determine which variables had the greatest effect on our resultant connectivity model. We adjusted resistance values assigned to 8 major variables included in our model: (1) natural landcover, (2) agricultural landcover, (3) rural landcover, (4) roads, (5) inland water, (6) elevation \ 200 m, (7) elevation between 200 and 1200 m, and (8) elevation [ 1200 m. We adjusted each resistance value assignment by both ± 20% and ± 50% in turn, while holding all other resistance values constant. We then reran the final connectivity model with the same buffer width and nodes. We compared these 32 resultant current maps against the baseline current map with Pearson's r to determine which variables, once adjusted, resulted in the greatest deviation from our baseline current map.

Uncertainty analysis
We conducted an uncertainty analysis to assess the degree to which changes in resistance values could potentially impact conservation decisions based upon our model. We quantified this impact by calculating an ''agreement score'' between scenarios-defined here as the degree to which a model identified the same 25 municipalities as having the highest mean current as other competing models. We calculated this agreement score for each scenario of adjusted resistance with the following methodology. We calculated mean current values across all municipalities included within our study area under each of the 32 scenarios tested in the sensitivity analysis; we treated the 20% and 50% sensitivity models as separate groups. We identified those 25 municipalities with the highest mean current under each scenario. Then, for each topranked municipality within each scenario, we calculated the percent of other scenarios in the group (including the original model) which also identified this municipality as having one of the highest 25 mean currents. We then averaged the score for each municipality within each scenario for the agreement score. Thus, if one sensitivity scenario has an agreement score of 80%, this indicates that, on average, the top 25 municipalities identified by this scenario were included in 80% of the other scenarios' top municipalities.

Characterizing Andean bear connectivity in Colombia
We performed several post-hoc analyses to characterize and better understand those areas of Colombia which were particularly conducive to connectivity across the Colombian Andes. First, we reclassified the cumulative current map into 5 classes using Jenks natural breaks (Jenks 1967). We categorized anything in the highest two classes as ''highly connective habitat'' and eliminated fragments that were less than 100 km 2 . We analyzed the proportion of this intact habitat that fell within different categories of protected areas according to a shapefile available from Colombia's National Natural Park Service (Parques Nacionales Nautrales de Colombia 2018). We also calculated the percentage of this habitat that either contained or was adjacent to agriculture as per the IDEAM 2014 landcover dataset.
Our sensitivity analysis revealed only slight changes to the final connectivity model following adjustments to resistance values. For example, a 50% increase in the resistance value for agricultural landcover (from 100 to 150) resulted in a connectivity model that was 99.521% correlated with the original model. Changes of ± 20% and ± 50% to resistance values for agricultural landcovers and the mid-elevational zone (between 200 and 1200 m) resulted in the greatest discrepancies with our final connectivity model (Fig. 5). In all scenarios, lower resistance values for each variable resulted in larger differences in final connectivity models than higher resistance values.
Models with the 50% adjustments to resistance values had, on average, 92.1% agreement among the 25 municipalities with the highest mean currents. The model with agricultural resistance reduced by 50% resulted in the least agreement; municipalities ranked by this model as having the 25 highest mean currents had only an 81.5% chance of appearing in another model's top 25 (Fig. 6b). The scenario with midelevation resistance reduced by 50% had more agreement with the others (88.25% agreement regarding top 25 municipalities). Uncertainty analyses with the 20% adjustments demonstrated an even higher degree of agreement (96.2% overall) (Fig. 6a).
We identified 41,487 km 2 as ''highly connective habitat'' (HCH) (this HCH largely corresponds to those areas depicted in blue in Fig. 3). Just over onethird of this HCH was contained within protected areas (38.5%). National natural parks represented the majority of this protected habitat (62.4%), but much Fig. 3 Cumulative current density map for Andean bears in Colombia derived from Circuitscape software and calculated with 78 nodes (3003 unique node pairs) randomly located at the perimeter of a 45-km buffer Fig. 4 a Mean cumulative current was significantly higher within those municipalities with records of Andean bear (Tremarctos ornatus) sightings (n = 155) than all other municipalities within the study area for which we found no records (n = 775); error bars represent 95% confidence intervals.
b Mean current density was significantly higher within 1 km of geographic coordinates of Andean bear (Tremarctos ornatus) locations (n = 235) than mean current density around random locations (n = 1175); error bars represent 95% confidence intervals of this area also fell within regional natural parks (17.3%), flora and fauna sanctuaries (10.6%), and integrated management districts (9.7%). Less than 1% of this protected HCH was contained with civil society reserves. Agricultural areas intersected 7.6% of the HCH. We identified the fifteen municipalities that contained the most km 2 of agriculture intersecting HCH (see Online Resource 3). The intersection of agriculture with HCH occurred most frequently in the department of Norte de Santander, but Boyacá also had several municipalities listed among those 15 municipalities with the highest amount of HCH in agriculture.

Discussion
Our model is the first approximation of landscape connectivity for Andean bears in Colombia, as well as the first connectivity model of the ecological integrity of the Colombian Andes. Our model suggests that landscape connectivity for the species is rather limited between the three cordilleras; this makes sense given the degree of human development in the intermountain valleys and lower elevations in these areas. Fortunately, a recent genetic analysis by Ruiz-García et al. (2020) demonstrated that despite these challenges for movement, Andean bear populations among the three cordilleras have maintained some reproductive and genetic connectivity (though genetic dissimilarities among the three cordilleras were detected). Our model suggests that close to 8% of the area contributing most to connectivity across the Colombian Andes coincides with agriculture such as pastures for cattle grazing or various possible crops. These agricultural areas may represent significant risk to Andean bears as long as rural Colombians view Andean bears as threats to their livelihoods. Our model allows us to identify those municipalities where such human-bear conflicts are likely to occur, and thus where strategies to mitigate the conflict are most urgently needed.
Though our model is coarse (1 km resolution), the outcomes from our model are at a scale that is relevant to Colombian conservation practitioners; many practitioners work at the scale of the municipality for project implementation. For example, environmental education and awareness raising is one of the four overarching objectives described in the National Program for the Conservation of the Andean bear in Colombia (Mayr Maldonado 2001). Some practitioners are pursuing this objective by hosting Andean bear festivals which draw in citizens from the entire municipality in which they are hosted (pers. observ.). Furthermore, as part of Colombia's transition to decentralized environmental governance, the regional environmental authorities were mandated by the Colombian constitution to coordinate and advise on environmental decisions made by municipalities within their jurisdictions (Blackman et al. 2004). Fig. 6 Degree of agreement among scenarios regarding the 25 municipalities with the highest mean current density after resistance scores were adjusted by 20% (a) or 50% (b); please note the different y axes. A score of 95% indicates that, on average, the municipalities ranked in that scenario as having one of the 25 highest mean current densities appeared in 95% of the other models' top 25 However, most of these authorities operate with extremely limited financial resources and personnel, meaning they must prioritize the extent to which they engage with different municipalities. Our model can help guide the selection of priority municipalities with which to build or improve these institutional relationships as part of their efforts to ensure adequate protection of key Andean bear habitat. Our model could also be used to determine where limited resources for fine-scale studies might best be invested. Uncertainty analyses demonstrated a relatively high degree of fidelity in the identification of those municipalities which contributed most to landscape connectivity. However, both the sensitivity and uncertainty analyses suggest that more information about how Andean bears move through agricultural landscapes would be helpful in improving our understanding of connectivity for this species in Colombia.
We used publicly available web records to demonstrate that Andean bear sightings occurred in municipalities identified as having higher current by our Circuitscape model. Precise location data from SiB and GBIF provided further support for our connectivity model. Overall, the validation technique appeared to work well: we found enough records to demonstrate model validity with statistical significance. The Andean bear is considered a charismatic species in Colombia (e.g., its image can be observed in many murals throughout Bogotá and rural Colombia), and perhaps because of this charisma, sightings of this elusive species often make local news. We chose not to use these data to inform model construction at the outset because there are some obvious, problematic biases. For example, news of Andean bears observed in national parks (where they are expected to occur) are usually not reported in media outlets, but Andean bears likely spend a large percentage of their time in these areas. Our objective was to create a countrywide connectivity model for a focal species, the Andean bear. Though the literature on landscape ecology has started to deemphasize single species approaches, many practitioners still focus their efforts on umbrella or flagship species for which there is strong public support and often more sponsors. Furthermore, the approach demonstrated here could easily be adapted to construct regional connectivity models of ecological integrity that could then be similarly validated with publicly available data and web records of charismatic fauna occupying these regions, even if these species were not the original focus of the study.
Data scarcity has long been recognized as an issue for landscape ecologists; necessity has led researchers to identify creative solutions for filling these gaps. For example, Newbold (2010) discussed the possibilities and limitations of using historical museum records for species distribution modeling. Calls for researchers to take advantage of citizen science data have become increasingly common (e.g., Devictor et al. 2010;Theobald et al. 2015;Kobori et al. 2016;Brown and Williams 2019). We add to these options the possibility for using local news and social media to gather coarse location information on charismatic fauna for the purposes of validating ecological models. These data are likely more abundant than we realize (as demonstrated by our success finding 186 place references through such means), and largely remain an untapped resource even in regions characterized by scarce ecological data and limited financial support for research. Indeed, this method may work even for species not typically considered ''charismatic''; for example, Balestrieri et al (2019) found over 200 publicly accessible records of marten (genus: Martes) locations in Italy and were able to use these data to guide the development of niche models.
We had no difficulties in applying the approach demonstrated by Koen et al. (2014) at our larger spatial extent (our study area was 417,895 km 2 ; their study area was 11,225 km 2 ). The updated Circuitscape in Julia ) allowed for fast computing across this large landscape in relatively short time (approximately 1 h for a model with 3003 unique node pairs). We found interesting discrepancies in the width and number of nodes required to address node-placement bias compared to the Koen et al. (2014) study. They suggested a buffer of at least 20% of the width of the study area would be required to remove node-placement bias. Our study area was approximately 500 km at its widest point; 20% would have been a buffer of at least 100 km. However, we found that a 45-km buffer was sufficient to mitigate node-placement bias. Additionally, we found more nodes were required in our study than they found were necessary for theirs (78 versus 18 nodes). This suggests there is more to this relationship yet to be understood, and we echo their suggestion for researchers to conduct their own sensitivity analyses with regards to requisite buffer width and number of nodes. There may be a mathematical relationship between the requisite number of nodes, buffer width, and measures of study area shape and size.