How to deal with ground truthing affected by human‐induced habitat change?: Identifying high‐quality habitats for the Critically Endangered Red Siskin

Abstract Species distribution models (SDM) can be valuable for identifying key habitats for conservation management of threatened taxa, but anthropogenic habitat change can undermine SDM accuracy. We used data for the Red Siskin (Spinus cucullatus), a critically endangered bird and ground truthing to examine anthropogenic habitat change as a source of SDM inaccuracy. We aimed to estimate: (1) the Red Siskin's historic distribution in Venezuela; (2) the portion of this historic distribution lost to vegetation degradation; and (3) the location of key habitats or areas with both, a high probability of historic occurrence and a low probability of vegetation degradation. We ground‐truthed 191 locations and used expert opinion as well as landscape characteristics to classify species' habitat suitability as excellent, good, acceptable, or poor. We fit a Random Forest model (RF) and Enhanced Vegetation Index (EVI) time series to evaluate the accuracy and precision of the expert categorization of habitat suitability. We estimated the probability of historic occurrence by fitting a MaxLike model using 88 presence records (1960–2013) and data on forest cover and aridity index. Of the entire study area, 23% (20,696 km2) had a historic probability of Red Siskin occurrence over 0.743. Furthermore, 85% of ground‐truthed locations had substantial reductions in mean EVI, resulting in key habitats totaling just 976 km2, in small blocks in the western and central regions. Decline in Area of Occupancy over 15 years was between 40% and 95%, corresponding to an extinction risk category between Vulnerable and Critically Endangered. Relating key habitats with other landscape features revealed significant risks and opportunities for proposed conservation interventions, including the fact that ongoing vegetation degradation could limit the establishment of reintroduced populations in eastern areas, while the conservation of remaining key habitats on private lands could be improved with biodiversity‐friendly agri‐ and silviculture programs.


| INTRODUCTION
One of the most promising applications of species distribution modeling (SDM) for conservation management is ranking areas by estimated habitat quality (Kramer-Schadt, Revilla, & Wiegand, 2005). This use of SDM assumes that areas with high probabilities of occurrence predict high-quality habitats (Franklin, 2010). However, species are not always present where high occurrence probabilities are predicted (Peterson et al., 2011). This mismatch between modeled predictions and field observations may result from problems with the SDM itself, such as conceptual errors (e.g., when models do not include biogeographical barriers or biotic interactions), limitations in variable selection (due a poor understanding of factors driving species distribution or use of outdated presence information with respect to environmental predictors used (Peterson et al., 2011). However, in other cases, this mismatch may be driven by anthropogenic processes such as increased poaching and overexploitation (Sánchez-Mercado et al., 2014) or land transformation due urban or agricultural development. For these reasons, field validation of SDM predictions is recommended; however, it is often not performed (Greaves, Mathieu, & Seddon, 2006).
The most frequent strategy for validating SDM predictions in the field is searching for the species of concern in areas with high values of predicted occurrence probabilities: Detections are interpreted as a confirmation of high-quality habitat (Bosso, Rebelo, Garonna, & Russo, 2013;Rebelo & Jones, 2010). However, if land cover transformation is gradual then, species detection is still possible where habitat has been partly degraded but not lost. Such "snapshot detections" of species occurrence may generate a misleading picture of relative habitat quality, which in turn could have disastrous consequences if, for example, the model is used to identify areas for the reintroduction of captivebred endangered species (Lahoz-Monfort, Guillera-Arroita, & Wintle, 2014). In such situations, a more nuanced, nonbinary approach to field validation is essential. On the other hand, lack of detection can be noninformative if the species is temporarily absent from high-quality habitat due to seasonal movements or has suffered strong declines due to habitat-unrelated threats such as poaching or other forms of wildlife extraction.
Clearly, validating SDM based only on species detections could be inappropriate in areas threatened by vegetation degradation and wildlife extraction. We therefore developed a new approach for field validation when a mismatch between model results and simple detection is likely. In this study, we used historical presence records for the Red Siskin (Spinus cucullatus), a low-abundance, Critically Endangered, and little-studied bird, to examine the combined utility of species distribution models and ground truthing via this new approach, to improve identification of key habitats for conservation interventions. The Red Siskin, a small Neotropical finch, is a particularly appropriate system in which to develop these alternative methods because it has been largely extirpated from its historic range across Venezuela, eastern Colombia, and Trinidad (Rodríguez, García-Rawlins, & Rojas-Suárez, 2015). Currently, the species persists in a few isolated populations within Venezuela and in a recently discovered small disjunct population in Guyana (Robbins, Braun, & Finch, 2003). The Red Siskin is listed as Endangered globally and Critically Endangered in Venezuela as a result of historic overexploitation for the specialized pet trade, and captive breeding and reintroduction have been recommended as management interventions; however, habitat loss is thought to be an important threat, although data are scan (Rodríguez-Clark et al., 2015).The natural habitat of the Red Siskin in Venezuela includes primarily tropical premontane humid and dry forests (Coats & Phelps, 1985) the latter of which are among the most threatened ecosystems globally and are endangered in Venezuela due to conversion for urbanization and agriculture (Rodríguez, Rojas-Suárez, & Giraldo Herández, 2010 (3) Where are remaining key habitats, or areas with both high historic occurrence probability and low landscape transformation? In addition, to demonstrating the usefulness of our approach for identifying key habitats for a threatened, elusive, and poorly studied species, we also aimed to examine implications for the threat status of the species and consider the consequences of our results for the design of effective strategies for reintroduction, including habitat conservation, which may be needed to achieve self-sustaining populations of the Red Siskin.

| Study area
Although the precise historical range of the Red Siskin is unknown, expert opinion can be used to identify an appropriate study area that is likely to contain this range. Experts agree that the distribution of this species in Venezuela is shaped principally by three factors: elevation, forest cover, and humidity (Coats & Phelps, 1985;Rivero Mendoza, 1983). The Red Siskin is thought to use a variety of habitats including dry deciduous woodland, mixed deciduous forest, evergreen forest, and the savanna-forest ecotone with daily and seasonal movements covering several kilometers from humid premontane forest into drier semideciduous forest and grassy clearings. Red Siskins use elevations from 200 to 1500 m (Coats & Phelps, 1985;Hilty, 2003;Rivero Mendoza, 2004). Based on these expert habitat descriptions, our approach for generating the study area polygon was based on estimating the Extent of Occurrence using the overlap of three environmental variables as proxies for factors describing suitable Red Siskin habitat: (1) elevation model (ELEV, 1 km of resolution; CGIAR Institute, 2010); (2) proportion of tree cover (TREE, 0.5 km of resolution; Hansen et al., 2002); and (3) aridity index (AI, 1 km; Zomer, Trabucco, Bossio, van Straaten, & Verchot, 2008).
Elevation values in Venezuela range from 0 to 4,382 m; we created a binomial layer defining as suitable (1) range from 200 to 1500 m, and values outside this range as unsuitable (0). Tree cover values in northern Venezuela ranged from 0% to 84%; because Red Siskins can use several vegetation types with different coverage, we used a broad threshold so that we could include a wider range of vegetation cover.
Thus, we defined as suitable tree cover values (TREE) from 10% to 100%. Aridity index (AI) represented precipitation availability over atmospheric water demand, the ratio of mean annual precipitation, and potential evapotranspiration, between the years 1950 and 2000. AI values ranged from 1,847 to 26,349, with higher values representing more humid conditions. We defined AI values from 2,000 to 6,500 to be suitable, which corresponded to semiarid and dry subhumid conditions (Trabucco & Zomer, 2009;Zomer et al., 2008).
For the study area definition and for the analysis of historical probability of occurrence (see below), we decided to change the spatial resolution of all predictive variables to the lowest spatial resolution (1 km) to minimize error due to spatial uncertainty in georeferencing historical distribution records (Wieczorek & Hijmans, 2004). For that, we aggregated (decreasing the resolution) using nearest-neighbor resampling method (by centering). This method is implemented in the function resample from raster package in R (Hijmans & Van Etten, 2012).
Finally, we overlapped all binomial layers to define a study area within which minimal habitat conditions for Red Siskins were met (90,060 km 2 ; Figure 1a).

| Historical probability of occurrence
We estimated the historical probability of occurrence (Ψ H ) for the Red Siskin within the study area using a compilation of historic and recent presence records and models based on maximum likelihood. We first made an extensive search for historic and current records of species presence from four sources: (1) Global Biodiversity Information Facility (GBIF) (2) national and international museums, (3) interviews with local ornithologists, and (4) a literature review. For all sources, we used keywords in English and Spanish related to all common and scientific names of the species (considering synonyms, alternative spellings, and subspecies) using the list compiled by Encyclopedia of Life Finally, we did a systematic review of the scientific literature using ISI Web of Knowledge and Google Scholar, using keys words related to the species and found seven published works (Coats, 1982;Coats & Phelps, 1985;Collar et al., 1992;López, 1991;Phelps & Phelps, 1963;Rivero Mendoza, 1983, 2004) from which we retrieved 332 records from 1867 to 1992. In total, we compiled 491 records of species presence from 1847 to 2013 ( Figure 1b).
We considered a record to be any discrete observation of one or more birds with a unique combination of the following information: (1) source, (2) coordinates, (3) observation year, (4) sex and development stage (adult or juvenile), and (5) quantity reported (number of individuals). For records without specific geographic coordinates, we used location descriptions (place names, geographic features, etc.) to assign latitude and longitude based on gazetteers (GIS Data Depo, DIVA GIS).
When contrasting coordinates were provided by each gazetteer, we calculated the mean value and error of latitude and longitude. If the error was larger than the cell resolution used to project our predictions (1 km 2 ) or if original coordinates had rounded two decimals, we discarded the record (Figure 1b).
To estimate the historical probability of occurrence (Ψ H ) for the Red Siskin, we used a maximum likelihood approach based on logistic regression to fit a species distribution model as a function of covariates, as implemented in R ("MaxLike"; Royle, Chandler, Yackulic, & Nichols, 2012). In addition to the three covariates used to define the study area, we also considered the 19 climatic variables in the WorldClim dataset as predictors (resolution 1 km; Hijmans, Cameron, Parra, Jones, & Jarvis, 2005). We evaluated redundancy and collinearity between all covariables using a hierarchical cluster analysis based on Pearson correlation. We defined a cluster as the group of variables with correlation <0.6 and selected one covariable for each cluster (Sarle, 1990). Thus, our complete model included the following eight variables: aridity index (AI), forest cover (TREE), mean diurnal temperature range (BIO02), isothermality (BIO03), annual temperature range (BIO07), mean temperature of the warmest quarter (BIO10), precipitation in the driest quarter (BIO17), and precipitation in the coldest quarter (BIO19). We applied a square-transformation to variables with considerable skew (TREE, BIO03, BIO10, BIO17, BIO18, and BIO19) and standardized all variables to a zero mean and unit variance, as recommended for the algorithm implemented (Royle et al., 2012).
To fit the occurrence model, we used only the 88 georeferenced presence records from 1960 or later, because our covariates were built with data from this date or later. The MaxLike approach assumed that detectability was constant over the study area (Royle et al., 2012).
Sampling effort was extensive enough over time (56 years) to have detected the species in the areas considered if it was present, leaving detection probability sufficiently uniform to meet this assumption.
However, MaxLike also assumed that sampling was random, which clearly was not the case, even though our dataset included almost all known sources of records. Thus, to surmount this problem, we applied a random sampling to the reports and repeatedly partitioned them into two independent subsets: an occurrence probability calibration subset and an occurrence validation subset (see details below; Franklin, 2010). The calibration subset was used to fit SDMs as described below and consisted of 66 reports (75% of the data). The validation subset (22 reports) was used to validate SDM performance. We repeated this two-way partitioning five times, which created replicates allowing us to directly evaluate data heterogeneity (Peterson et al., 2011).
To select the "best" MaxLike occurrence probability model, we then fit different combinations of the eight covariates described above to each of the five replicate calibration data subsets. Our first model (mdl1) contained linear terms for all eight variables. The second (mdl2) included only climatic variables, and the third model (mdl3) was the most reduced, including only forest cover and aridity index, which have been proposed by experts to be the most important variables affecting Red Siskin occurrence (Rivero Mendoza, 2004). The "best" model was considered to be the one that both converged and had the lowest AICc in the most replicates (Table 1; Burnham & Anderson, 2002).
We used evaluate and threshold functions from dismo package in R (Fielding & Bell, 1997) to (1) select the replicates of the best model with the best performance to built spatial prediction, and (2) select the threshold of historic probability of occurrence at which Red Siskin presence is the highest (T ΨH ). Model prediction was evaluated based on correlation coefficient (cor), Area Under the Receiver Operator Curve (AUC), and maximizing the sum of sensitivity and specificity (maxSSS) using as pseudoabsences 88 points that were randomly selected from the northern part of the country, but outside of the study area for the Red Siskin. We selected the replicates 1, 2, and 5, which had the highest values of AUC, cor, and maxSSS to built the spatial prediction (Table 2). We used the mean value of statistic "max kappa" (predicted value at which kappa is highest; Liu, Berry, Dawson, & Pearson, 2005) of selected replicates as criteria to set T ΨH (0.743).
To built the spatial prediction, we used the predict function of maxlike package (Royle et al., 2012) and a raster stack of the same predictive variables disaggregated at a resolution of 250 m to produce a map that were comparable with our other predictions (see below).

| Ground-truthing and vegetation time series analysis
We randomly selected 90 points within the study area ( Figure 1c) to perform ground truthing of habitat suitability. At each point, we walked transects of 1.5 km, which we laid on the roads nearest to the point. Along each transect, we stopped every 500 m, resulting in 270 "evaluation locations" in total (90 random points * 3 stops per point). However, several of these evaluation locations were in areas that were either inaccessible or in areas with high risk to the personal security of field teams. We, therefore, systematically discarded 79 locations present in dangerous or inaccessible areas, resulting in a total of 191 evaluation locations. on geographic characteristics (elevation, slope), and vegetation type (forest, shrubs, and pasture) that are considered important according to expert accounts and own field experience (Coats & Phelps, 1985;Rivero Mendoza, 1983, 2004. Locations with "excellent" suitability were those with mosaic of forests surrounded by shrubs, located at medium elevations (500-800 m) and steep slopes (>60%). Locations T A B L E 1 Statistical support (AICc values), and convergence status for three models of red siskin occurrence, fit with MaxLike to five replicate data subsets  In order to relate the subjective habitat suitability classification with the measured EVI phenology, we fit a random forest classification model (RF). We implemented RF in the randomForest package in R (Liaw & Wiener, 2002 Therefore, we applied a matrix of ordered weights to recalculate the OOB (Piccarreta, 2008).
We visualized the spatial distribution of habitat with current optimal suitability using the predict function of randomForest package EVI also allows the meaningful comparison of seasonal and interannual changes in vegetation growth and activity. We therefore used the time series of EVI data to describe changes in the mean and variance of EVI values within the period studied at each evaluation location. Changes in local vegetation did not occur simultaneously for all evaluation locations. Thus, we used the method proposed by Chen and Gupta (2000)

| Identifying key habitats
Finally, we overlapped the historic probability of occurrence (from our SDM) and current habitat suitability predictions (from our RF) to identify "key habitats," or areas with both high historic occurrence probability and excellent habitat suitability. To identify key habitats, we multiplied the values of both predictions to generate an overlap index that ranged from 0 to 2.248, with 0 indicating low historic probability of occurrence and low habitat quality and values close to 2 indicating high values for both conditions. We also performed an extinction risk assessment using the criterion of population reduction based on an estimated decline in Area of Occupancy (AOO) and habitat quality (criterion A2c; IUCN, 2012). We estimated AOO based on the number of cells containing key habitats and calculated the percentage of reduction in this area with respect to the historically suitable area for Red Siskins (i.e., areas with Ψ H > 0.743). We also calculated AOO when key habitats were defined to include both "good" and "excellent" habitat suitability, to take into account uncertainty in Red Siskin habitat use.

| Spatial distribution of habitat with current optimal conditions
The majority of the evaluation locations (85%) had a substantial reduction in mean EVI values over the last 15 years (Figure 2). These changes were similar across different categories of habitat quality (chi-square = 4.707; df = 3; p = .195). For 60% of the evaluation locations, this change occurred before 2011. Locations classified with optimal habitat suitability ("good" and "excellent") consistently had historic mean EVI values >0.4. Although for these same locations, current mean EVI values were substantially lower, they were still above 0.4. Locations classified with suboptimal habitat suitability (acceptable and poor) were more heterogeneous, with EVI values from 0.3 to 0.6, but mostly below the mean values of the optimal habitats ( Figure 2).
The overall corrected classification error rate of the RF model was 33.2%, with lowest classification error for "excellent" (22.1%) largest for "acceptable" (57.1%). Habitats classified as "good" were predicted in a wide area within the study area (35,494 km 2 , 39% of study area), while "excellent" habitat (3,127 km 2 , 3% of study area) were clustered in the western part of the country, in the lowlands of the Sierra de Perijá and along the southern slope of the Cordillera de Mérida ( Figure 3a). Other small and more dispersed blocks of "excellent" habitat were predicted in the center-west as well as in the east. Habitat with suboptimal conditions ("poor" and "acceptable") was focused in three large blocks (55,307 km 2 , 61%) in the west, center, and east of the country.

| Historical probability of occurrence
Presence records before 1960 (25 records) were located mostly in the western part of the country, while more recent sightings (94) were evenly distributed across the center and west (Figure 1b).
The best model for historic probability of occurrence was mdl3 (containing only the aridity index and forest cover). Alternative models containing other climatic variables did not converge due to scarcity of records (Table 1). Our estimates of historical occurrence probability were based on the three replicates of mdl3 that had good predictive accuracy ( Table 2). The area with the highest occurrence probabilities (Ψ H > 0.743) covered 20,696 km 2 and was concentrated toward the center and north of the Coastal Cordillera. In the west, small fragments with high probability were observed toward the south, along the north slope of the Andean Cordillera and Serranía de Portuguesa, and along the eastern coast of Lake Maracaibo (Serranía del Empalado). To the east, there were also discontinuous fragments around the Turimiquire Massif and west of the Araya Península (Figure 3b).
Of the area with the best historically suitable area for Red Siskins (Ψ H > 0.743), only 4,686 km 2 (23%) was protected in a Venezuelan national park. The most valuable unprotected habitats were in the western region, including areas in the northern Sierra de Perijá, along the west coast of Lake Maracaibo and in the mountains of Falcón and Lara. In the northeast, there was also a wide continuous area with high probabilities on unprotected lands (Figure 3b).

| Key habitats
Key habitats, defined as areas with both high historic occurrence probabilities and currently "excellent" suitability, covered just 976 km 2 and occurred in the western and central regions (northern end of the Sierra de Perijá, lowlands of Sierra de El Empalado, and the Coastal Cordillera), forming small blocks (Figure 3c). Only three small If key habitats were expanded to include both "good" and "excellent" areas, their area increased by an order of magnitude, to 12,274 km 2 .  , 2012). Given the scarcity of records throughout this range, it is furthermore possible that the actual area occupied by the species is far less than the area available. This is likely because we hypothesize that in addition to habitat loss, we suspect the Red Siskin has suffered what is known as a "high-abundance-biased" or HAB decline (Rodríguez, 2002): Individuals were likely removed from the geographic range not randomly or evenly from across the range, but rather in a way that was biased toward high-abundance areas. This is because trappers seem to have been specialized, and interested in this species in particular, and so likely searched for and trapped it precisely in the areas that they were most likely to find it. Our analysis provides the first quantitative evidence that in addition to overexploitation, land transformation may also be driving the extirpation of Red Siskins in Venezuela, and also reveals that ongoing habitat transformation could limit the establishment of reintroduced populations there (Figure 3a and c). The distribution of key habitats corresponded well to the conservation status described by Coats and Phelps (1985) Although habitats with optimal conditions have thus far retained some forest cover (EVI values above 0.4), the ongoing degradation observed suggests that remaining blocks with optimal habitat conditions could also be degraded in the short-to-medium term, reducing the availability of suitable habitat for this critically endangered species even more. Parks; Figure 1a and 3c). Interestingly, the objective delineation of key habitats described here also helps bring into focus opportunities for conservation action. For example, between Henri Pittier and Macarao National Parks, a potential corridor includes 200 hectares currently covered with shade coffee farms. This agroforestry habitat currently faces an uncertain future due government price restrictions that make standard coffee production unprofitable (SUNDE, 2015). The Red Siskin Initiative has proposed to apply a proven market-based approach, Bird Friendly Coffee ® certification (BFC) to these shade coffee farms, which would qualify their products as a specialty coffee, free of price restrictions (Philpott, Bichier, Rice, & Greenberg, 2007). BFC certification could be a means to protect and improve the shade coffee farm habitat present in this corridor, preserving a potential reintroduction site for Red Siskins that is also prime habitat for migratory birds.

| Model accuracy
The Random Forest model is used here to transform a subjective evaluation of habitat quality into an spatial index of habitat suitability for conservation planing. This application assumes that the categories suggested by experts are indeed predictive of the occurrence and viability of the species in the field, and the selected variables are good indicators of the expert ranking. This is, however, a difficult task, given the inherent uncertainty in expert opinions and the natural variability in environmental conditions. The accuracy of the RF model in predicting habitat suitability categories based on EVI time series suggested a moderate overall performance. However, the model was better at discerning optimal than suboptimal habitat conditions. The largest classification errors occurred in habitats with "acceptable" categories, which covered a wide type of vegetation conditions. These errors could be due to a lack of understanding of habitat requirements for such a rare and little-studied species, which could generate an underestimation of the amount of habitat available to Red Siskins. However, this error could also reflect the capacity of Red Siskins to use transformed habitats, such as ecotones of dry deciduous woodlands, shrubby grasslands, and pastures (Robbins et al., 2003). where agricultural production and other human activities are segregated from areas managed for biodiversity conservation. Clearly, achieving integration between human activities and conservation objectives for Red Siskins in Venezuela requires a more detailed understanding of temporal and spatial patterns of species habitat use. Even more importantly, this will require building capacity for rural communities to adopt biodiversity-friendly land management (e.g., water source protection, healthy soil management, sustainable agroecology) and the promotion of policies that encourage them (Brussaard et al., 2010).
Our approach has proven useful for identifying key habitats for a threatened and poorly sampled species and also to monitor temporal and spatial trends in vegetation transformation within these habitats.
In the case of the Red Siskin, this approach is only the first step toward identifying suitable habitat for reintroduction, which should be refined with additional research focused on breeding and feeding ecology, seasonal movements, and the spatial distribution of poaching risk (e.g., Sánchez-Mercado, Asmussen, Rodríguez-Clark, Rodríguez, & Jedrzejewski, 2016).

Funds for this research were provided by the Smithsonian Grand
Challenges Consortia, the Instituto Venezolano de Investigaciones Científicas (IVIC), and Birdfair/RSPB. We are grateful to M. Lentino, D. Ascanio, C. Sharpe, J.G. León, G. Rodríguez, F. Escola, and L.
Moran for providing their observational records of Red Siskins.