Introduction

The Harpy Eagle (Harpia harpyja) has a broad distribution from southern Mexico to northern Argentina yet has been extirpated across much of its former range, particularly Mexico and much of Central America (Ferguson-Lees and Christie 2001, Vargas González et al. 2006). Due to recognized declines, its global conservation status is Near Threatened (International Union for Conservation of Nature 2017), while its national conservation status in many countries is Endangered (e.g., Panama: Ministry of the Environment of Panama [MiAMBIENTE] 2016; Belize: Meerman and Clabaugh 2010).

Behaviors that render Harpy Eagles largely inconspicuous (e.g., individuals rarely soar, still-hunt from perches inside the canopy, and are often silent; Stiles and Skutch 1989) and low breeding densities (Piana 2007, Muñiz-López 2008, Vargas González and Vargas 2011) make them difficult to locate, thus presenting challenges for observing breeding biology and behavior. Pairs build nests in canopy-emergent trees (Schulenberg 2020), where they raise a single young that requires 2.5 to 3 yr to reach independence (Vargas González and Vargas 2011). Breeding occurs at low elevations (e.g., 300–500 masl) in primary and mature secondary forests (Vargas González and Vargas 2011, Schulenberg 2020), a habitat that provides suitable nesting sites and abundant large, arboreal prey needed for reproduction (Brown 2001, Geise et al. 2004, Please et al. 2016). Harpy Eagles are often considered to be sensitive to anthropogenic disturbances including habitat alteration (Thiollay 2007, Aguiar-Silva et al. 2012) and human persecution (e.g., shooting due to fear or curiosity; Trinca et al. 2008, Curti and Valdez 2009, Muñiz-López 2017), both of which are factors that potentially limit breeding near human population centers (Vargas González 2008). At the micro-habitat scale, Harpy Eagles show specific preferences in habitat selection at nesting sites (Giudice et al. 2007, Piana 2007, Vargas González et al. 2014) that probably derive from the distribution and availability of prey and high perches, and offer ease of mobility between forest strata to move and hunt, features that contribute to breeding success and survival of the offspring during the dependency period (Vargas González 2008).

Although we possess basic knowledge about Harpy Eagle breeding biology, habitat selection at landscape levels remains unquantified. Given the increasing pressures on Harpy Eagle habitat and conflicts with humans, the development of a breeding habitat distribution model could help in conservation planning and prioritization of conservation actions.

Species distribution models relate environmental variables to species occurrences to gain an understanding of ecological drivers or to help predict habitat suitability across landscapes (Elith and Leathwick 2009, Kramer-Schadt et al. 2013). Numerous studies have documented the practical usefulness of predictive models for habitat suitability (Sara 2014, Cardador et al. 2015, Hamilton et al. 2015, Maleki et al. 2016, Munro et al. 2016, Gallardo and Vilella 2017) for the conservation and management of species. Some authors have indicated that the use of these models for threatened species is essential because they aid in prioritizing areas for conservation and facilitate the development of management recommendations (Moisen et al. 2006, Seoane et al. 2006).

Our aim was to develop a predictive model for Harpy Eagle breeding habitat as a tool to aid in the species' conservation. We included in the analysis natural and anthropogenic variables taken from 27 nesting sites in the western Pacific region of Darien Province, Republic of Panama, and used MaxEnt to identify and quantify suitable breeding habitat. We then compared proportions of suitable habitat among watersheds and landscapes with varying levels of administrative oversight and protection.

Study Area and Methods

Study Area. Our study area is located in Darien Province, Panama, between 7°33.0′ and 8°18.0′N latitude and 77°24.6′ and 78°12′W longitude, totaling 5636 km² (Fig. 1). The area experiences marked dry (January–April) and wet (May–December) seasons, with an annual average rainfall of 1372 mm and average temperature of 28°C (Instituto Nacional de Estadística and Censo de Panamá 2016).

Figure 1.

Map delineating the study area, forested habitat, and nest sites used in the habitat distribution model for Harpy Eagle habitat suitability in Darien Province, Panama.

Darien is considered part of the Chocó biogeographic region, which extends from southern Panama to northern Ecuador and is globally recognized as one of the world's most biologically diverse (Gentry 1986, Dinerstein et al. 1995). Three of the five climatological regions of Panama (classified according to Köppen 1936) are found in Darien, which creates niches for biodiversity, and consequently makes Darien an ecological region of international conservation importance (Myers et al. 2000, Olson and Dinerstein 2002). A large part of Panama's biodiversity is found in Darien ecosystems, with 54% of Panama's vertebrate species and 47% of endemic species occurring there (Tosi 1971, United Nations Development Program – Ministry of Economics and Finance of Panama 2003). Four protected areas occur in Darien: Darien National Park, Chepigana Forest Reserve, Serranía de Bagre Biological Corridor, and Punta Patiño Wetland.

Nesting Sites. We follow Franke (2017) in defining nesting sites as the specific location on the landscape where nests (the actual structure where eggs are laid) are located. We georeferenced (10-m accuracy) 41 Harpy Eagle nesting sites monitored as part of our long-term (2000–2017) research, primarily in the northern section of the study area (Fig. 1). Due to access issues, local permit acquisition, and safety issues, the northern area of the study region was sampled more heavily than the southern. This information was collected as part of a long-term program for Harpy Eagle conservation and research conducted by The Peregrine Fund in Panama and supported by local people, who offered anecdotal information that we systematically verified to confirm Harpy Eagle nesting sites. We verified whether sites were occupied or unoccupied through direct observation, and used as evidence of breeding the observation of an incubating adult or a nestling, fledgling, or dependent juvenile < 1 yr old (identified by plumage coloration). We followed the recommendation of previous authors and excluded alternate nesting sites from analyses to reduce the chance of pseudoreplication (Álvarez-Cordero 1996, Vargas González and Vargas 2011). In Panama, the use of alternate nests is uncommon and usually results from a mitigating factor (e.g., habitat alteration near the original nest), and we expect that excluding alternate nests from analysis is biologically and statistically justifiable.

Nesting-site Variables. Harpy Eagle home ranges in Panama average 16–24 km² (Vargas González and Vargas 2011). Based on prior research, we considered five natural and anthropogenic nesting site-level variables within 250 m (the resolution of our coarsest spatial data) of the nest and assumed to be important in selection of nesting habitat by Harpy Eagles (Vargas González 2008, Vargas González and Vargas 2011, Vargas González et al. 2014): slope, elevation, distance to rivers, distance to large communities (>500 inhabitants), and percent forest cover within 250m. We also included two additional forest cover variables that represent selection of primary breeding habitat at larger landscape scales (within 1 km and within 4 km of the nest), for a total of seven variables (Table 1). All variables of interest and their derived products were projected in EPSG: 32617, WGS 84, with a resolution of 250 m.

Table 1.

Mean, standard deviation, and range of the seven variables used in the predictive habitat suitability analysis for Harpy Eagle nest locations, background points, and the entire study extent.

We used Program R (R Core Team 2016) to process the spatial data. To extract elevation and calculate slope, we used a digital terrain model (DEM) from the Shuttle Radar Topography Mission project (Shuttle Radar Topography Mission 2000, US Geological Survey 2006). Large communities were delineated according to the criteria of the Panama Statistics and Census Institute. Rivers were digitized (maximum scale 1:20,000) from parameters obtained at the Ministry of Public Works and the Cartographic Institute of Panama. The land cover layer was obtained through Panama's Ministry of the Environment (MiAMBIENTE 2003). We then aggregated primary and mature secondary forested land cover into one forested type excluding all non-forest land cover types (water, agriculture, flooded vegetation, urban, and mangrove). Lastly, to calculate percent of forested habitat within 1 km and 4 km of nests, we used a moving window method (raster package in R; Hijmans 2016). This method retains the 250-m resolution; however, each cell or “window” represents 1 km or 4 km surrounding each focal cell. Eagles might actually be selecting breeding sites at even larger scales than 4 km, however, we were limited to a distance of 4 km from the nest because our land cover data did not extend into Colombia, which borders our study region along the southeast (Fig. 1).

Modeling Habitat Suitability. We modeled relative probability of occurrence, hereafter referred to as habitat suitability, using known nesting sites and available sites (i.e., background data) as a function of seven nest- and landscape-level covariates. We used MaxEnt 3.4.1 (Phillips et al. 2006, Phillips and Dudík 2008) in the R dismo environment (Hijmans et al. 2017). MaxEnt is a general-purpose machine learning method with a predictive performance consistently competitive with the highest performing presence-only methods (Elith et al. 2006, 2011), especially for small sample sizes (Hernandez et al. 2006).

A key assumption of presence-only models is that sampling is either random or representative throughout the landscape (Yackulic et al. 2013). However, presence-only data typically result from studies where presence records are spatially biased toward better surveyed areas (Kramer-Schadt et al. 2013). Since spatial bias generally results in environmental bias, the difference between occurrence records and background sampling may lead to inaccurate models (Phillips et al. 2009, Kramer-Schadt et al. 2013, Merow et al. 2013, Yackulic et al. 2013). Ideally, sampling bias would be incorporated into the model using a “bias file” representing the relative sampling effort (Kramer-Schadt et al. 2013) or “target-group” sampling, using presence locations of taxonomically similar species to estimate effort, under the assumption that those surveys would have recorded the focal species had it occurred (Phillips and Dudík 2008, Merow et al. 2013). However, because nests were most often found opportunistically by local villagers who reside in the northern region of Darien, no survey effort or additional species occurrences were collected. All nests found by villagers were then validated by researchers to confirm species identification and eagle breeding status. To account for the strong bias in survey effort in the northern region, we chose background points to match the sampling bias. We randomly selected 10,000 background points from a buffer of 5 km surrounding the nesting sites. Means of the seven variables used in our analyses were relatively similar among nests, background points, and the entire study area (Table 1). Knowing we were extrapolating to unsampled areas, we created multivariate environmental similarity surfaces (MESS; Elith et al. 2010, 2011) to visually depict a measure of similarity between the training sample and the study extent, therefore, indicating where extrapolation will occur (Fig. 2). As expected, dissimilarity decreased as the extent of availability points got larger. However, dissimilarity was low even for the 5-km buffer, with only 17% of the cells indicating a negative dissimilarity value and only 3% showing high (dark blue) negative values, meaning increased extrapolation in those areas (Fig. 2).

Figure 2.

Multivariate environmental similarity surface (MESS) depicting the similarity between three training samples (5-km buffer, 10-km buffer, and the study extent) for Harpy Eagles in Darien Province, Panama. Red indicates high similarity (positive values) between samples and extent, while blue (negative) shows the areas with the highest differences (i.e., areas where extrapolation is likely to occur).

Additionally, we spatially thinned the 41 nests by a distance of 1.24 km, the mean distance reported for alternative nests (Vargas González and Vargas 2011). This was done to reduce individual bias due to pooling nests across years and some pairs (33%) having multiple nests (Vargas González and Vargas 2011). Spatially thinning reduced the sample size to 27 nesting sites and only two nesting sites were from the same pair; however the nests were found in different years (Fig. 1).

To generate models with maximum predictive ability, we considered models that exhausted potential covariate combinations with restrictions on model complexity. We only considered additive combinations of linear and quadratic covariate relationships because of the small sample size. All seven covariates were included, regardless of multicollinearity. However, only the three forest cover variables were significantly correlated using Pearson's correlation coefficient (>0.7: Dormann et al. 2013) with each other (  Table S1 (349_rapt-54-04-11_s01.pdf)). Machine learning algorithms can receive correlated predictors and many studies have shown little difference in results when correlated variables are excluded (Elith et al. 2011, Kramer-Schadt et al. 2013). However, care should be taken when interpreting response curves of correlated predictor variables (Boria 2014). To reduce overfitting, MaxEnt allows a user-defined parameter to adjust the strength of the regularization penalty that penalizes the model proportional to the magnitude of coefficients, in turn shrinking many coefficients to zero, essentially removing them from the model (Phillips et al. 2006, Merow et al. 2013). Models were tuned across a range of complexity by varying the regularization multiplier from 0.5 to 4.0 in stepwise increments of 0.5. This resulted in a candidate set of eight models. We used the package ENMeval (version 0.3.0; Muscarella et al. 2014) to tune the models.

ENMeval also allows various methods to partition occurrence locations into testing and training groups for k-fold cross-validation. We chose the “block” method, which partitions the data according to latitude and longitude lines, and divides the data into four folds of equal numbers. Models are then run iteratively using k – 1 bins for training and the remaining one for testing (Muscarella et al. 2014). We chose this method because it accounts for spatial autocorrelation between testing and training data resulting from biased sampling and has shown greater transferability across time and space than other methods (Wenger and Olden 2012, Muscarella et al. 2014).

Model Evaluation. The final model was selected using two evaluation metrics: omission rates and area under the curve (AUC). Omission rates calculated using the 10% calibration omission threshold (OR) reflect overfitting to the calibration data when over 10% test omission (Shcheglovitova and Anderson 2013, Boria 2014). AUC of the receiver evaluates the model's ability to discriminate conditions at withheld occurrence locations from those at the background samples (Muscarella et al. 2014, Radosavljevic and Anderson 2014). An AUC = 0.5 indicates discrimination that is no better than random and AUC = 1 indicates prefect discrimination (Fielding and Bell 1997, Pearce and Ferrier 2000). Models with AUC > 0.7 have good discriminatory power (Hosmer and Lemeshow 1989, Kramer-Schadt et al. 2013). However, Yackulic et al. (2013) states AUC for presence-only modeling is a relative value for comparing performance between different models created from the same data and should not be subject to a singular threshold. Therefore, we chose the model with the lowest omission rate and highest AUC. Additionally, we tested the predictive ability of our final model using an out-of-sample method. We binned our predicted habitat suitability values into six bins using a quantile method (Morris et al. 2016). This method divides the continuous values of 0–1 into six bins of equal number of data points (i.e., cells). We chose six bins for two reasons; first, ten bins are excessive with only 14 out-of-sample nests, and second, we wanted to combine an equal number of bins to create a low-med-high suitability gradient for our conservation areas analysis. We then calculated a Spearman-rank correlation coefficient between frequency of the excluded data points within each bin and the bin rank (1 = low to 6 = high). A strong positive correlation indicates the model has good predictive accuracy (Boyce et al. 2002).

Prioritization of Areas for Conservation. We obtained protected area boundaries for the Republic of Panama from MiAMBIENTE and overlaid these and the study regions' watersheds to locate and quantify areas with suitable habitat. We examined habitat suitability at the watershed scale because MiAMBIENTE is currently using watersheds as a method to systematize conservation actions. For both protected areas and watersheds analysis, we reduced the six habitat suitability bins to three, resulting in a low–medium–high suitability gradient.

Results

The 27 Harpy Eagle nests used in the analysis were found in areas of high forest cover at the level of nesting site and surrounding landscape (Table 1). Nesting sites had a mean elevation of 138 masl, mean slope of 10°, and on average were found 1.3 km from a river and within 16.9 km of a large community (Table 1). We selected the model with the highest AUC and lowest OR from the candidate model set of eight models that varied in regularization coefficients. The regularization coefficient of the top model was MaxEnt's default of 1. The model fit the training data and discriminated nesting sites from background locations relatively well with AUC values of 0.76 and 0.74, respectively. The model had an OR of 0.18, indicating very slight overfitting (i.e., >0.10). For comparison, models fit with regularization coefficients of 0.5–2.5 had the same OR values, and coefficients of 3–4 had values >0.22. We transformed the raw model output using a complementary log-log (cloglog) transform (Phillips et al. 2017; Fig. 3).

Figure 3.

Predicted binned cloglog habitat suitability (0 = low to 1 = high) for Harpy Eagle nesting sites in the Darien Province, Panama, study area.

Our out-of-sample predictive performance test indicated the final model showed high predictive performance. The Spearman-rank correlation between bin rank (1–6) and frequency of out-of-sample nests was strongly positive (r = 0.95, P = 0.005).

The final model contained all predictor variables except forest cover within 250 m of the nesting site (Table 2). The most influential variables were forest cover within 1 km of the nesting site, elevation, and distance to rivers (Table 2). Habitat suitability values were nearly identical at forest cover (1 km) from 0–40%, whereas the increase was nearly exponential from 40–100% (Fig. 4). Contrary to forest cover 1 km, habitat suitability decreased abruptly in areas where elevation was greater than 300 masl (Fig. 4). Distance to rivers showed a positive relationship with habitat suitability. However, habitat suitability was never <0.5 (Fig. 4). Slope, forest cover at 4 km, and distance to large communities added little to the model (Table 2). Habitat suitability increased with increasing slope; the highest suitability was between approximately 10–50° (Fig. 4), which corresponded to the range of slopes in the study extent (Table 1). Forest cover at 4 km and distance to large communities both had a slightly negative relationship with habitat suitability, and although the relationship was negative, habitat suitability never fell below 0.5, as also found for distance to rivers (Fig. 4).

Table 2.

Percent contribution of each predictor variable to the final model of habitat suitability for Harpy Eagles in Panama.

Figure 4.

Partial dependency plots showing the marginal effect of the predictor variables and the predicted habitat suitability for Harpy Eagles, Darien Province, Panama. Percent forest at the nest site (forest cover 250 m) is not shown here because it did not contribute to the model.

By collapsing six bins into three, we created a habitat suitability gradient from low to high (Fig. 5). Of the 5636 km² of land area analyzed, 1888 km² (34%) showed highly suitable breeding habitat, 1868 km² (33%) ranked medium suitability, and 1880 km² (33%) ranked low habitat suitability for the Harpy Eagle according to the model (Table 3). Of the 14 out-of-sample nest sites, 12 were located in highly suitable habitat, two in medium, and none in low suitability.

Figure 5.

High – medium – low habitat suitability gradient for Harpy Eagles overlaid with watershed boundaries (labeled), and four protected areas in Darien Province, Panama. Cells were classified using a quantile method. Watersheds are identified by the names of their respective rivers (Tucutí, Sambú, and Juradó), or by rivers that mark their boundaries (Tucutí and Sambú, Sambú and Juradó), a naming protocol used by Panama's Ministry of the Environment (MiAMBIENTE).

Table 3.

Estimates of the extent of breeding habitat for the Harpy Eagle in five watersheds of Darien, Panama. Breeding habitats are classified as high suitability, medium suitability, or low suitability, based on our modeling, and estimates of area and proportion of suitable habitat for each watershed are shown. Watersheds are identified by the names of their respective rivers (Sambú, Tucutí, and Juradó), or by rivers that mark their boundaries (Tucutí and Sambú, Sambú and Juradó), a naming protocol used by Panama's Ministry of the Environment (MiAMBIENTE).

Of the five watersheds we examined separately, the Tucutí watershed was the largest with 1838 km² (33% of total land area), followed by the Sambú, the watershed between the Sambú and Jurado rivers, the watershed between Tucutí and Sambú rivers, and the Juradó watershed with areas of 1496 km² (26%), 1070 km² (19%), 891 km² (16%), and 341 km² (6%), respectively (Table 3, Fig. 5). The basin with the greatest area of highly suitable habitat was Tucutí (665 km², 35% of total) followed by Sambú (462 km², 25% of total), the watershed between the Tucutí and Sambú rivers (399 km², 21%), the watershed between Sambú and Juradó rivers (329 km², 17%), and Juradó (33 km², 2% of total; Table 3, footnote 1). When considering both high and medium suitable habitat together, the pattern remained the same: Tucutí watershed with the most suitable habitat (36% of total), followed by Sambú (26%), Tucutí and Sambú (21%), Sambú and Juradó (15%), and Juradó (2%). Of four protected areas located in the study area, 833 km² (15% of total) had highly suitable habitat (Fig. 5). The protected area with the greatest proportion of highly suitable habitat was Darien National Park with 21% (392 km²), followed by the Serranía de Bagre Biological Corridor with 12% (232 km²), and finally, the Chepigana Forest Reserve with 11% (209 km²; Table 4, footnote 1). The same pattern holds true when considering both high and medium suitable habitat together; Darien National Park contains the most (29% of total), then the Serranía de Bagre Biological Corridor (11%), and lastly the Chepigana Forest Reserve (9%). The fourth, Punta Patiño Wetland, located within the northern region of the La Palma watershed, contained <1% suitable habitat and is not further addressed here.

Table 4.

Estimates of the extent of breeding habitat for the Harpy Eagle in three protected areas of Darien Province, Panama. Breeding habitats are classified as high suitability, medium suitability, or low suitability, based on our modeling. Estimates of area and proportion of suitable habitat for each protected area are shown.

Discussion

The observed specificity of suitable habitat for nesting sites suggests that the Harpy Eagle has a reduced ecological amplitude. In particular, we observed a positive association with forest cover in the landscape (1 km), a greater likelihood of occurrence below 300 masl, and distances of <1 km to rivers. Habitat specificity observed in this study, combined with the eagles' nearly clumped distribution at landscape and regional levels, their low abundance, and the progressive destruction of forests in the study region, argues for recognition of the Harpy Eagle as a conservation-dependent species (Global Raptor Information Network 2017), and designation of their breeding habitat as vulnerable and requiring sustainable management.

We developed the first predictive habitat distribution model for Harpy Eagle breeding habitat. We found the model had good predictive ability, and therefore can be used to guide conservation actions benefitting Harpy Eagles. Variables included in our habitat distribution model suggest that Harpy Eagle nesting sites occurred in areas with high forest cover. Forest cover within 1 km of nesting sites made the greatest contribution to the models, whereas forest cover within 4 km had a minor influence, and forest cover within 250 m had no influence. We speculate that high forest cover within 1 km of nesting sites may provide complex forest structure and sufficient prey densities needed by juvenile eagles to improve flight and hunting skills during their nearly 2-yr post-fledging dependency period. Limited data on juvenile dispersal indicate that most movements in the first 2 yr after fledging occur within approximately 1 km of the nest (Muniz-Lopez 2012). The dependency on forests for nesting is consistent with results from prior research in Panama (Vargas González and Vargas 2011), Peru (Giudice et al. 2007), and Ecuador (Muñiz-Lopez 2008).

Elevation and slope are two important physio-graphic factors that influence suitability of Harpy Eagle nesting sites. Although elevations in the study area ranged from 0 to 5922 masl, nesting sites were found at low elevations, predominantly approximately 100 masl, and elevation was the second-most important variable in the model. The biological reasons Harpy Eagles nest at lower elevation are less clear. Preferred prey, especially sloths (Bradypus variegatus and Choloepus hoffmanni) and kinkajous (Potus flavus), occur from low to middle elevations (Muñiz-Lopéz 2008, Vargas González and Vargas 2011), and diversity and abundance of Harpy Eagle prey may be greater at lower elevations (Pires et al. 2000). Additionally, forest structure at lower elevations may be more suitable for breeding and hunting. A common trend in Neotropical forest structure is a decrease in forest height with elevation, with most canopy volume found in mid-levels of lowland forests (Asner et al. 2014, Girardin et al. 2014). Lower elevations have tall trees suitable for nesting; in addition, lower elevations provide a forest profile of emergent canopy above the dense mid-story, a structure suitable to the perch-and-ambush hunting method favored by Harpy Eagles. Despite published records of Harpy Eagles observed at high elevations (Álvarez del Toro 1980), our observation of nesting at low elevations is consistent with findings from Peru (Giudice 2007), Guyana (Rettig 1978), and Ecuador (Muñiz-Lopez 2008). Slope was only minimally influential in informing the model, and we observed Harpy Eagle nests on slopes averaging 10°. Perhaps the association between elevation and slope influences microclimatic factors (e.g., temperature, light intensity, and humidity) favorable to incubation, hatching, or nestling-rearing.

We found that Harpy Eagle nests were more likely to be closer to rivers (mean distance 1.25 km) compared to available locations. Harpy Eagles use open spaces around rivers to observe, chase, and catch prey (Vargas González 2008). Forest openings associated with rivers also promote the growth of plant species (e.g., Cecropia spp.) with leaves and fruits favored by primary (kinkajous and sloths) and secondary (primates) prey species (Franco-Rosselli and Berg 1997, Kays 1999, Ramirez et al. 2011). Interestingly, the model indicated a positive relationship, meaning suitable habitat increases with increasing distance to rivers, which may have been an artifact of the small range of distance values due to the extensive river system in the study area.

Given our approach, and the conservative distance from nests that available points were selected, distance to large communities (mean 17 km) was not an important predictor in our model (Table 2). Prior research in Panama and Venezuela (Álvarez-Cordero 1996) and in Peru (Piana 2007) documented some tolerance to anthropogenic pressure in their breeding habitat, but the level and type of disturbance that breeding Harpy Eagles can tolerate remain unquantified. We suspect that probability of incidental shooting increases directly with the number of humans living near nesting sites (Trinca et al. 2008), and that increasing distance from nesting sites to population centers provides a buffer against persecution. Further, habitat disturbance increases with size of settlements, and diminishes with distance from settlements, such that intact forests needed for nesting are likely to occur at greater distances from large settlements than small ones. As human populations grow and spread in Darien, the availability of intact forests free from persecution will increasingly become a limiting factor for Harpy Eagle reproduction.

Eighty-six percent (n=12) of the nesting sites used in the out-of-sample validation process were found in highly suitable habitat, while the remaining two validation nests were found in medium suitable habitat and none in low suitability. Based on area requirements and density estimates (Vargas González and Vargas 2011), the highly suitable breeding habitat in the study region can harbor a population of between 76 and 113 breeding Harpy Eagle pairs. We recognize that additional improvements to the model could be derived by including other data that condition the selection of nesting sites, such as the abundance and diversity of prey and human hunting pressure, data that are unavailable at the present time. We envision practical uses for the model, such as testing hypotheses regarding productivity at sites with high and low habitat suitability.

Despite the ability of our model to discriminate high quality nesting habitat, we still observed one Harpy Eagle pair (of the 27 used to create the model) nesting in low suitability habitat and three in medium suitability habitat. We suspect that the existence of such nesting sites at the time of analysis may have been due to fidelity to ancestral nesting sites (Viter 2013) despite recent anthropogenic changes to surrounding habitat. Since 2000, Panama has experienced significant deforestation (Food and Agriculture Organization of the United Nations 2011). To explore forest loss within the landscape around the nests used in the model building, we calculated the percent forest loss within 1 km of nests using the Hansen et al. (2013) global forest loss layer. We chose the 1-km scale because it proved the most important scale in the model (Table 2), with the highest habitat suitability between 40% and 80% forest cover. We found that nest sites experienced highly variable forest cover loss, with a median loss of 3.4% (range = 0.1–45.5%) between 2000 and 2018; our monitoring in 2017 indicated 88% of the monitored nests were unused (i.e., alternative nests were not known, and we assume that breeding attempts did not occur at 88% of territories in 2017). Whether the abundance of unused sites was due to forest loss or another human pressure is unknown. As stated above, data on human persecution would be invaluable in determining nesting-site selection. We also note that converting the continuous model output to a gradient of low to high may not be biologically meaningful for the Harpy Eagle.

In terms of protected areas, Darien National Park contained 392 km² of highly suitable habitat, approximately equal to the combined areas of Chepigana Forest Reserve (209 km²) and Serranía de Bagre Biological Corridor (232 km²), two administrative units that share a common border (i.e., geographically and biologically they are a single forest). This observation highlights the fact that lands outside Darien National Park in the northern portion of the study area are equally important for conservation of Harpy Eagles. This is a concern for Harpy Eagle conservation in Darien, because extractive forestry and other land use practices that damage forest cover and, according to our model, will reduce habitat suitability, are permitted under authority of MiAMBIENTE (Panamá 1960, 1980, 1994, 1995a, 1995b) in both Chepigana and Serranía de Bagre. There are initiatives for forest use and exploitation in these regions, for which rapid ecological assessments have been developed (González et al. 2011, Fuenmayor 2011) that do not describe the importance of these habitats for the species; a situation that potentially puts the largest known population of Harpy Eagles in Central America at risk (Vargas González et al. 2006).

The Harpy Eagle is the National Bird of Panama (Panamá 1995, 2002) where its conservation status is classified as Endangered (MiAMBIENTE 2016). Considering the spatial distribution of suitable breeding habitat for Harpy Eagles in and outside protected areas with varying levels of forest exploitation, improvements to current initiatives in place for conservation are warranted. We recommend reevaluation of the management categories of the Chepigana Forest Reserve and the Serranía de Bagre Biological Corridor, given the importance of their ecosystems for the reproduction and stability of Harpy Eagle populations.

Implications for Conservation. Our predictive model clearly identified areas of importance for the conservation of Harpy Eagles. The greatest conservation implication of this model lies in its practical utility to identify the landscape where a population that is in danger and vulnerable is concentrated due to the accelerated and dynamic transformation of the habitat. Our model allows managers to prioritize targeted conservation initiatives, such as environmental education, community development programs, and sustainable community-based activities, to help curb the deterioration of natural ecosystems, and to combat Harpy Eagle persecution. Likewise, it allowed us to determine that the areas with high suitability lacked governmental legislation restricting the use of forest resources, a situation with negative implications for the eagle and the underlying biodiversity that coexists in the natural environment where it is found. Although the effect of the extirpation or population reduction of the Harpy Eagle in the ecosystems where it coexists with other species has not been evaluated, the importance of top predators to maintain the integrity and equilibrium in the ecosystem is well known (Bierregaard 1998, Beschta and Ripple 2011). Finally, this study offers a replicable methodology to be applied and tested in other regions.

Supplemental Materials (available online).  Table S1 (349_rapt-54-04-11_s01.pdf): Pearson's correlations coefficients for all predictor variables used in the model.