Using surrogate species and MaxEnt modeling to prioritize areas for conservation of a páramo bird community in a tropical high Andean biosphere reserve

ABSTRACT

The páramo grassland ecosystem in the high Andes requires definition of spatially explicit and large-scale priorities for bird conservation, especially for lands outside of protected areas. Using surrogate species such as habitat specialists (e.g., páramo specialists), endemic species, or threatened species to identify potential habitat can support conservation decisions. We used MaxEnt to generate habitat suitability models for eighteen surrogate bird species in a high Andean biosphere reserve, the Macizo del Cajas (MCB). The areas of estimated suitable habitat range from 115,500 to 312,700 ha; the proportion predicted inside the national system of protected areas of MCB ranged from 12 to 29 percent. The consolidated map (i.e., the coincidence of predicted pixels for all eighteen surrogate bird species) predicted an area of 94,800 ha, with 70 percent outside the national system of protected areas. The models reveal that there are large suitable areas of MCB potentially occupied by these species outside of protected areas, especially along the eastern flank of the study area. MaxEnt models as a proxy facilitate the recognition of locations of possible habitat suitability for páramo-specialized birds and therefore identify areas of the MCB outside of the protected areas that should be considered for more formal protection.

KEYWORDS:

• Ecuador
• Macizo del Cajas
• páramo ecosystem
• habitat suitability
• protected areas

Introduction

Ongoing human modification of ecosystems increasingly confines biodiversity to disturbed areas with reduced amounts of natural habitat (Foley 2005; Gardner et al. 2009). Human modification of natural areas is an important driver of biodiversity loss, particularly in regions with high levels of endemism and high concentrations of threatened species (Ceballos et al. 2015; Waldron et al. 2017). Identifying priorities for conservation efforts is crucial, particularly in biodiversity hotspots (Myers et al. 2000). The tropical Andes is a priority conservation area because it harbors not only one of the greatest concentrations of species restricted to a particular area or ecosystem but also one of the highest concentrations of threatened bird species in South America (Stotz et al. 1996; Stattersfield et al. 1998; Devenish et al. 2009), thus qualifying as a hotspot of bird diversity and endemism (Myers et al. 2000). Within the tropical Andes, the páramo grassland ecosystem of the northern high Andes occurs across Venezuela, Colombia, Ecuador, and northern Perú, with patchy distribution between Costa Rica and Panamá. These ecosystems form an ecologically isolated enclave with extreme daily temperature fluctuations throughout the year, mist cover, and high rates of solar irradiation (Hofstede and Llambí 2020). Organisms living here are adapted to extreme environmental conditions, resulting in important processes of speciation, endemism, and remarkable biodiversity (Hofstede and Llambí 2020); a mark of this distinctive heterogeneous ecosystem is its mosaic of habitats (Neill 1999). For instance, at a fine scale, the páramo grassland ecosystem is characterized by a combination of vegetation types such as native bunch grasses (tussock páramo grassland), taller semi-open shrubland (shrubby páramo), marshes with bog plants in more humid areas (cushion páramo), and numerous patches of Polylepis forest (Neill 1999; Sklenár and Ramsay 2001; see Figure 1 for example photos). These habitats harbor important regional bird diversity across the páramo ecosystem (Jiménez-Rivillas et al. 2018; Astudillo et al. 2020). Birds that exhibit a high affinity or exclusivity to these habitats are considered habitat specialists (because they occupy mainly one or a few of these Andean habitats; Fjeldså and Krabbe 1990; Kessler et al. 2001; Devenish et al. 2009; Jiménez-Rivillas et al. 2018).

Figure 1. Example photos of important habitats for birds across the páramo ecosystem in the Macizo del Cajas Biosphere Reserve, southern Andes of Ecuador. We present two examples per habitat (upper and lower panels): (a) tussock páramo grassland, (b) cushion páramo, (c) shrubby paramo, and (d) Polylepis woodlands.

The dependence of bird species to high-altitude ecosystems, such as páramo grasslands, makes them important indicators for conservation prioritization across the high Andes (Myers et al. 2000; Sala et al. 2000; Lees et al. 2020). As with grasslands worldwide, bird species adapted to páramo habitats have been negatively influenced by agricultural expansion and conversion of native grassland to pastures for livestock grazing (Sala et al. 2000; Foley 2005; Norment, Runge, and Morgan 2010; Newbold et al. 2016; Hannah et al. 2020). In the Andean region, at least 70 percent of high-altitude habitats have been strongly affected by human activities (Critical Ecosystem Partnership Fund 2021). The high-altitude Andes are also expected to suffer dramatic habitat modifications due to climate change, with greater seasonality considerably affecting habitat configurations with vegetation adaptable to more extreme climate regimes becoming more dominant (Foster 2001). This in turn will influence the availability of microhabitats for mountain habitat specialists (Foster 2001; Barnett, Adam, and Lettenmaier 2005), such as birds that occur across the páramo grassland ecosystem. The páramo grassland ecosystem is limited by the montane forest line (~3,400 m.a.s.l.) and the mountaintops (~4,500 m.a.s.l.) or the permanent snow line (~5,000 m.a.s.l.; Neill 1999). This is key, because the lower limit is predicted to rise due to climate change, leaving a narrower band of suitable environmental conditions for páramo habitats and their associated birds (Dybala, Gardali, and Eadie 2013; Campos-Cerqueira et al. 2017). Both land use and climate change lead to habitat disturbance that will affect the distribution and suitability of habitat for mountain birds.

In the tropical high Andes, spatial variation is strongly influenced by topographic and climatic patterns such as elevation and precipitation (Jiménez-Rivillas et al. 2018; Udy et al. 2021). This variation results in environmental heterogeneity (Udy et al. 2021). Heterogeneous environments offer more niche partitioning between species, resulting in high levels of biodiversity and coexistence among organisms (e.g., Ben‐Hur and Kadmon 2020). For instance, the complex topography contributes to both connectivity (e.g., corridors across massifs, continuous mountains ranges) and isolation (e.g., patches at the top of the mountains, deep valleys) of páramo habitats, resulting in a large number of different vegetation associations, with each one representing a particular community assembly (Hofstede and Llambí 2020). There is more remaining suitable habitat for fauna adapted to high-altitude habitats in areas with less accessible topography, which leads to less human activity (Gale 2000; Sylvester et al. 2017). Furthermore, the elevation gradients in the Andes create a high variation of microclimates with more seasonality in precipitation as elevation increases, and this pattern shows regional delineated precipitation over small spatial scales (e.g., Ballari et al. 2018), which can drive ecological divergence by promoting specific conditions for habitat-specialized birds (e.g., Benham and Witt 2016). All of these phenomena promote greater complexity of structure and composition of vegetation (i.e., habitat heterogeneity) with consequently greater bird diversity across the Andean region (Poulsen and Krabbe 1998; Willig and Presley 2016; Wallis et al. 2016; Sevillano-Ríos and Rodewald 2017; Astudillo et al. 2019). Because the presence of these habitats is driven by these topographical and climatic variations, understanding these patterns is likely to help identify areas of habitat suitability for Andean birds. Therefore, knowledge of the relationship between environmental variables (e.g., topography and climate) and avian diversity is crucial for identifying habitat suitability, particularly for habitat-specialized species such as páramo birds (Norambuena and Van Els 2021).

The maximum entropy algorithm (MaxEnt) has been widely applied to estimate habitat suitability based on environmental conditions and occurrence data to improve conservation efforts through estimating high habitat suitability for focal species (Hijmans and Graham 2006; Gallagher, Hughes, and Leishman 2013; Fastré et al. 2020). MaxEnt is an advanced modeling technique used to predict important areas via estimating environmental characteristics at presence sites (occurrence) and at background locations (study region); the locations may be selected randomly or in a regular pattern (Guillera-Arroita, Lahoz-Monfort, and Elith 2014). MaxEnt has been used in Ecuador to identify important biodiversity areas, including within the Macizo del Cajas Biosphere Reserve (southern Ecuador; Lessmann, Muñoz, and Bonaccorso 2014; Cuesta et al. 2017). However, in previous studies, páramo birds were lumped within a generalized perspective of multitaxa biodiversity conservation (e.g., Cuesta et al. 2017), with little attention applied to specific conservation concerns such as páramo-specialized species particularly in páramo regions with important levels of diversity and endemism such as Macizo del Cajas Biosphere Reserve (Jiménez-Rivillas et al. 2018). For instance, violet-throated metaltail (Metallura baroni) is an endangered and endemic hummingbird that only occurs within the highlands of Macizo del Cajas (Tinoco et al. 2009; Freile et al. 2019). Other birds, such as tit-like dacnis (Xenodacnis parina) and stout-billed cinclodes (Cinclodes excelsior), have larger local populations in the country as well as in the highlands of Macizo del Cajas (Astudillo, Tinoco, and Siddons 2015). In addition, within the Macizo del Cajas, the MaxEnt algorithm has only been applied for bird conservation for the potential distributions of two specific species, violet-throated metaltail (Tinoco et al. 2009) and Andean condor (Vultur gryphus; Astudillo et al. 2016). All of the aforementioned species are associated strongly with the páramo grassland ecosystem and are considered focal for conservation (Astudillo, Tinoco, and Siddons 2015). Here, these focal species could act as a surrogate for identifying habitat suitable for most páramo bird species. The concept of surrogate species is to use a small (i.e., subset) but representative number of species to contribute to solving conservation problems (Caro and O’Doherty 1999). However, no work to date has used MaxEnt modeling with surrogate species that are specifically important for the conservation of the páramo bird community as whole. Here, we used a subset of the community that presents relevant characteristics for conservation priorities of the entire páramo bird community as surrogate species (Caro and O’Doherty 1999; Caro 2010). We considered surrogate species of the páramo bird community as those characterized as habitat specialist (i.e., high affinity to páramo habitats), endemic species, and/or threatened species; identifying suitable habitat for these surrogate bird species may improve recognition of key conservation areas for páramo suitable for the entire avian community.

The Macizo del Cajas Biosphere Reserve is located in the southern high Andes of Ecuador and includes two core areas with protected status (i.e., Cajas National Park and Quimsacocha National Recreation Area). The páramos of Macizo del Cajas are classed as a biogeographic district with specific levels of endemism (Jiménez-Rivillas et al. 2018). However, 84 percent of the páramo landscape in the biosphere reserve is beyond the limits of these protected areas where livestock farming and burning to promote pasture is common (Astudillo et al. 2018; Barros et al. 2020). Páramo habitats are also under ongoing pressure from road infrastructure, expansion of croplands, and introduction of exotic tree plantations (Astudillo 2014; Aguilar et al. 2019). Thus, combining the overlapping habitat suitability estimations of surrogate bird species (by consolidating habitat suitability models) is an important tool to direct conservation priorities in the páramo grassland ecosystem throughout the Macizo del Cajas Biosphere Reserve, a conservation hotspot. To do so, we used four years of páramo bird surveys to model habitat suitability of endemic, threatened, and páramo-specialized bird species as surrogate species to identify key areas for bird conservation. We use this approach to identify páramo areas with potentially high concentrations of these surrogate bird species to reveal priority conservation areas, with particular relevance to lands currently lacking formal protection.

Methods

Study area

The Macizo del Cajas Biosphere Reserve (MCB) is mainly located in the provinces of Azuay and Cañar, southern Ecuador. The MCB is bounded to the north by the Cañar River and to the south by the Jubones River (Figures 2 and 3) and covers an area of 976,601 ha. The entire MCB contains several western Ecuadorian ecosystems, including coastal, humid and arid lowlands, the Andean west slope, and the high Andean cordillera. This study was conducted across the páramo ecosystem in the high Andean cordillera of MCB (> 3,400 m.a.s.l., 2°55′25″ S, 79°21′57″ W). The páramos of the study area are considered a particular biogeographic district within the páramo biogeographic province of the northern Andes (Jiménez-Rivillas et al. 2018) with high levels of diversity and endemism.

Figure 2. Habitat suitability maps of eighteen surrogate bird species using the MaxEnt algorithm. The first eighteen maps show the individual species predictions (see bird codes in Table 1) and the respective occurrences (dots). The two final maps show the species aggregation map (see methods) created by summing all reclassified individual distribution models and a map of the central coordinate of mountain transects across the páramo landscape (brownish polygon) of the Macizo del Cajas Biosphere Reserve (yellowish polygon), southern Andes of Ecuador. All distribution models were selected via tenth percentile or higher (see Table 1). The red polygons represent the national system of protected areas. The arrows on the first map indicate the western and eastern flanks of the study area (see study area description) omitted in further maps to improve interpretability.

Figure 3. Consolidated map of the area predicted to contain all eighteen surrogate bird species in the páramo landscape of Macizo del Cajas Biosphere Reserve (yellowish polygon), southern Andes of Ecuador. Area shown in green is where all predicted habitat suitability of individual surrogate species from MaxEnt models overlap (pixel value = 18; see Methods). The brownish polygon is the distribution of páramo grassland ecosystem (>3,400 m.a.s.l.). The arrows indicate the western and eastern flanks of the study area (see study area description).

Table 1. List of surrogate bird species of the Macizo del Cajas Biosphere Reserve, southern Andes of Ecuador.

Scientific name	English name	Species code	Habitat affinity	Region of endemism	Threat	AUC	The top four variable contributions (%)	Tenth percentile	Total predicted area (ha)	Proportion of predicted area outside of protected areas (%)	Number of occurrences, training records, and testing records
Oreotrochilus chimborazo	Ecuadorian hillstar	ORCH	P, Ps			0.997	Elevation = 56.4	0.27	222,500	83	47 (training = 33; testing = 10)
							Precipitation seasonality = 19.1
							Isothermality = 12.5
							Mean diurnal range = 4.6
Aglaeactis cupripennis	Shining sunbeam	AGCU	P, Ps			0.998	Elevation = 53.2	0.3	312,700	88	41 (training = 23; testing = 7)
							Precipitation seasonality = 15.4
							Isothermality = 15
							Mean diurnal range = 9.5
Metallura baroni	Violet-throated metaltail	MEBA	P, Ps, F	Central Andean páramo	EN	0.997	Elevation = 49.4	0.28	165,200	78	37 (training = 31; testing = 10)
							Precipitation seasonality = 17.8
							Isothermality = 11.6
							Mean temperature of wettest quarter = 7.3
Chalcostigma stanleyi	Blue-mantled thornbill	CHST	P, Ps		VU	0.997	Elevation = 58.6	0.29	130,800	71	129 (training = 45; testing = 15)
							Precipitation seasonality = 20.5
							Isothermality = 10.2
							Mean diurnal range = 3.4
Leptasthenura andicola	Andean tit-spinetail	LEAN	P, Ps			0.989	Elevation = 61.9	0.32	172,700	78	145 (training = 51; testing = 16)
							Precipitation seasonality = 18.2
							Isothermality = 10.1
							Mean diurnal range = 4
Asthenes flammulata	Many-striped canastero	ASFL	P, Ps			0.995	Elevation = 54.6	0.33	182,400	80	395 (training = 72; testing = 23)
							Precipitation seasonality = 21
							Isothermality = 9.8
							Mean temperature of wettest quarter = 9.3
Cinclodes excelsior	Stout-billed cinclodes	CIEX	P	Central Andean páramo		0.996	Elevation = 48.7	0.25	165,100	78	38 (training = 26; testing = 8)
							Precipitation seasonality = 19.1
							Mean temperature of wettest quarter = 14.2
							Isothermality = 13.2
Cinclodes albidiventris	Chestnut-winged cinclodes	CIFU	P			0.992	Elevation = 58.6	0.29	217,700	82	273 (training = 73; testing = 24)
							Precipitation seasonality = 20
							Isothermality = 10.1
							Mean temperature of wettest quarter = 4
Grallaria quitensis	Tawny antpitta	GRQU	P, Ps, F			0.995	Elevation = 53.4	0.35	216,300	82	338 (training = 79; testing = 26)
							Precipitation seasonality = 19.7
							Isothermality = 11.5
							Mean temperature of wettest quarter = 7.4
Agriornis montanus	Black-billed shrike-tyrant	AGMO	P			0.99	Elevation = 54.2	0.38	197,700	81	163 (training = 63; testing = 21)
							Precipitation seasonality = 19.2
							Isothermality = 11.1
							Mean temperature of wettest quarter = 8
Muscisaxicola alpinus	Plain-capped ground-tyrant	MUAL	P			0.997	Elevation = 47.4	0.26	193,700	80	50 (training = 36; testing = 11)
							Precipitation seasonality = 21.4
							Isothermality = 12.4
							Mean temperature of wettest quarter = 10.6
Ochthoeca fumicolor	Brown-backed chat-tyrant	OCFU	P, Ps, F			0.996	Elevation = 55.4	0.36	213,900	82	139 (training = 51; testing = 17)
							Precipitation seasonality = 18.4
							Isothermality = 11.9
							Mean diurnal range = 6.36
Cnemarchus erythropygius	Red-rumped bush-tyrant	CNER	P, Ps, F		NT	0.997	Elevation = 58.8	0.32	147,500	76	38 (training = 18; testing = 6)
							Precipitation seasonality = 15.4
							Isothermality = 10.5
							Mean diurnal range = 5.9
Cistothorus platensis	Grass wren	CIPL	P			0.995	Elevation = 51.7	0.37	200,200	71	260 (training = 63; testing = 20)
							Precipitation seasonality = 18.2
							Isothermality = 12.2
							Mean temperature of wettest quarter = 9.2
Xenodacnis parina	Tit-like dacnis	XEPA	Ps, F		EN	0.998	Elevation = 57.5	0.39	115,500	73	45 (training = 21; testing = 7)
							Precipitation seasonality = 15.8
							Isothermality = 9
							Slope = 7.1
Catamenia inornata	Plain-colored seedeater	CAIN	P			0.995	Elevation = 55.7	0.32	215,100	82	71 (training = 36; testing = 12)
							Precipitation seasonality = 17.1
							Isothermality = 12.7
							Mean diurnal range = 5.6
Catamenia homochroa	Paramo seedeater	CAHO	P			0.997	Elevation = 50.8	0.29	188,900	81	53 (training = 33; testing = 11)
							Precipitation seasonality = 17
							Mean temperature of wettest quarter = 11.7
							Isothermality = 11
Geospizopsis unicolor	Plumbeous sierra finch	PHUN	P			0.996	Elevation = 57.1	0.34	203,100	81	258 (training = 73; testing = 24)
							Precipitation seasonality = 19.1
							Isothermality = 10.2
							Mean temperature of wettest quarter = 5.5

This list includes habitat affinity, endemism criteria, as well as IUCN threat category for each species. We present the MaxEnt results for predicting their distribution including the AUC for the model, the four most important variables in each model, the value of the tenth percentile and the number of occurrences per bird species (see methods). We also present the total predicted area and the proportion of that area outside the national system of protected areas. Scientific and English names follow the South America Classification Committee (Remsen et al. 2021). Habitat affinity follows the classification suggested by Astudillo, Tinoco, and Siddons (2015) and Barros et al. (2020). The endemism classification follows Stattersfield et al. (1998). The threat categories follow the red list of birds of Ecuador (Freile et al. 2019). The habitat affinity codes are as follows: P = tussock grass páramo; Ps = shrubby páramo; F = Polylepis forest. The threat categories are as follows: NT = near threatened; VU = vulnerable; EN = endangered.

The MCB was designed around two core protected areas administrated within the national park system across the high Andean cordillera: Cajas National Park (PNC; 2°50′45″ S, 79°14′33″ W) and Quimsacocha National Recreation Area (ARQ; 3°00′45″ S, 79°14′12″ W). These protected areas represent 3 percent of the total area of MCB and contain 16 percent of the reserve’s páramo grassland ecosystem (Barros et al. 2020). In addition, the MCB considers adjoining areas around both protected areas as buffer zones (40 percent of the total surface of MCB; Barros et al. 2020). These buffer zones are a combination of public, private, and communal territories located across the páramo grassland landscape of MCB. In the study area, ~90 percent of the vegetation is páramo grassland (Neill 1999; Figures 1–3). This páramo grassland is an open habitat characterized by native tussock grass species Figure 1(a) in association with cushion páramo plants in more humid areas Figure 1(b), as well as shrubby páramo with native woody shrubs in semi-open habitats Figure 1(c); Astudillo et al. 2018). Forest fragments dominated by Polylepis tree species are found throughout the páramo grassland ecosystem Figure 1(d); Astudillo et al. 2020). Consequently, the regional páramo is a complex matrix of páramo habitats (Neill 1999) with greater heterogeneity (botanically and structurally) concentrated in valleys or on steep slopes, particularly in areas with difficult access for human activity (Astudillo et al. 2018; Barros et al. 2020). However, there is also an evident natural homogenization of páramo habitat with less woody plant cover and increasing proportion of tussock grass cover as elevation increases, especially at elevations ≥3,900 m (Astudillo et al. 2018, 2019). The average annual rainfall ranges from 1,200 to 1,500 mm and generally shows a bimodal rainfall with the greatest rainfall occurring between March and May but with a second, less intense peak in rainfall from September to February; the driest period is from June to August (Celleri et al. 2007; Ballari et al. 2018). However, the eastern flank of MCB (Figures 2 and 3) is primarily influenced by the eastern Amazon rainfall regime, which mainly follows a unimodal regime (Campozano et al. 2016; Ballari et al. 2018), which results in more tussock grass cover in association with humid cushion bogs and fewer woody plants Figures (a,b); Barros et al. 2020). In contrast, the western flank of MCB (Figures 2 and 3) facing the Pacific lowlands mostly follows the bimodal rainfall regime, with relatively lower mean annual precipitation (Campozano et al. 2016; Ballari et al. 2018). This western climate is associated with more native woody plant cover, fewer humid areas, and less tussock grass cover as well as cushion bogs Figure 1(c); Barros et al. 2020). The average monthly temperature ranges from 5°C to 12°C on both flanks.

Species occurrence and MaxEnt modeling

We used information generated by four years of bird surveys performed between February 2016 and November 2019 within MCB. Surveys consisted of 106 strip transects (each 1 km in length) randomly installed across the páramo grassland landscape of MCB (Figure 2). Each transect was separated from each other by at least 500 m in open and semi-open habitats (e.g., tussock páramo grassland, shrubby páramo, cushion páramo) at elevations from 3,450 to 4,150 m. Bird counts along the transects were always conducted 15 minutes after sunrise until 11 o’clock. Walking at constant speed (~0.5 km h⁻¹), we recorded all birds seen and heard within 50 m of each side of the observer, excluding flyovers. Two observers performed each transect. The same observer always conducted the counts, while the second observer noted exact GPS coordinates for each bird recorded for the species occurrences. Each transect was sampled at least three times a year to encompass the climatic seasonality of the study area. Although we monitored the entire bird community, we selected only eighteen surrogate species to include in the MaxEnt models.

Surrogate species are a select subset of a community whose attributes are of conservation concern. In this study, our surrogate species were chosen to help us to identify areas of páramo beyond the limits of the protected areas of MCB that these birds may occupy and therefore contribute to more effective bird conservation at the community level. Thus, we included páramo bird species that represent conservation concerns (i.e., endemic and threatened). Endemism was determined based on the list of endemic bird species by Stattersfield et al. (1998), and threatened species status was based on the red list of bird species of Ecuador (Freile et al. 2019). Habitat affinity has also been shown to be an effective tool for identifying conservation priorities throughout the study area (e.g., Latta et al. 2011; Astudillo et al. 2018, 2020; Barros et al. 2020). Therefore, we also incorporated bird species that are resident in the study area and are specialized in or use the páramo habitats (i.e., have an affinity for tussock páramo grassland, shrubby páramo, cushion páramo, and/or Polylepis woodland habitat) as surrogate species (Table 1). Hence, the groups of habitat affinity are páramo, birds that prefer open areas associated with tussock grass as well as cushion paramo; shrubby páramo, birds that prefer semi-open areas with taller vegetation characterized by native woody shrubs; and Polylepis forest, birds that prefer woodlands dominated by Polylepis plant species (Astudillo, Tinoco, and Siddons 2015; Barros et al. 2020). We assigned a specific GPS coordinate for each bird observed, and the occurrence data for the eighteen surrogate species ranged from 37 to 395 georeferenced observations (Table 1).

We applied the maximum entropy algorithm (MaxEnt v3.4.0; Phillips and Dudík 2008) to obtain a picture of environmental characteristics at presence sites and at background locations. We selected MaxEnt over other algorithms because it has been tested extensively and has been found to perform suitably as an advanced modeling technique for predicting relative suitability regardless of how the background sample is specified (Merow, Smith, and Silander 2013; Guillera-Arroita, Lahoz-Monfort, and Elith 2014). MaxEnt calculates the ratio between the probability distribution describing characteristics of pixels at which species occur and the probability distribution describing the environmental characteristic as a whole (Phillips and Dudík 2008; Merow, Smith, and Silander 2013). The approach of maximum entropy to estimate the suitable habitat of species is equivalent to minimizing the environmental space (Phillips and Dudík 2008; Merow, Smith, and Silander 2013; Morales, Fernández, and Baca-González 2017). Therefore, MaxEnt does not estimate occurrence probabilities; rather, it represents relative suitability according to environmental data (Guillera-Arroita, Lahoz-Monfort, and Elith 2014). Because our objective is finding potential areas for priority conservation of páramo birds within the MCB and our occurrence data are derived from field observations across páramo habitats that cover large areas of the páramo region of MCB (Figure 2), we considered presence-only data as the most appropriate approach (i.e., yielding a ratio of environmental data where the species are observed with the environmental data in the background).

We modeled the eighteen surrogate species occurrences in the same run using climatic and topographic information for the whole of Ecuador as environmental variables. The environmental information has a resolution of 30 arcsec (~1 km²). We used eight environmental variables, including precipitation amount of warmest quarter (kg m⁻²), coefficient of variation of precipitation (precipitation seasonality), precipitation amount of wettest quarter (kg m⁻²), mean temperature of diurnal range (°C), mean temperature of wettest quarter (°C), isothermality (relation between the diurnal temperature range and annual temperature range), elevation (m.a.s.l.), and slope (%). Climatic information was downloaded from CHELSA-Free climate data (Karger et al. 2017, 2020). The CHELSA climatological data incorporates corrections in relation to orography (e.g., valley exposition, wind fields, boundary layer height), and therefore the climate information, particularly the data associated with precipitation, is more accurate (Karger et al. 2017, 2020; Beck et al. 2020) and is preferable for mountain systems such as the high Andes. Topographic variables were generated from a digital terrain model based on contour lines (1:3,000 scale). We included the same type of environmental data used previously for páramo specialized birds (i.e., Andean condor, violet-throated metaltail) because they represented a relatively high contribution to previous MaxEnt outputs in the study region (Tinoco et al. 2009; Astudillo et al. 2016). Though these two species represent widely different natural histories and ecological roles, similar environmental variables were found to be important and thus we assume that these same variables will provide useful information at the scale of páramo-specialized birds.

We evaluated the correlations among environmental variables (Appendix). Potential problems for highly correlated variables were tested by evaluating their contribution (Table 1) on the MaxEnt solution. For the most important variables for MaxEnt, only elevation is correlated with mean temperature of wettest quarter (r = −0.99), mean temperature of diurnal range (r = 0.52), and slope (r = 0.66). In addition, we used the regularized training gain data via jackknife checkbox as a second method to evaluate the contribution of the variables. The jackknife compares the contribution of each variable in isolation and their contribution in conjunction with the rest of the variables (Phillips and Dudík 2008). The models showed that mean temperature of the wettest quarter, mean temperature of diurnal range, and slope (correlated variables) have low contributions to the MaxEnt solution in both percentage contribution and the regularized training gain (via jackknife checkbox). The rest of the correlated variables have very little contribution to the MaxEnt solution (via percentage value of contribution and jackknife checkbox). Elevation showed a relatively higher percentage contribution (Table 1) as well as higher regularized training gain via jackknife (both in isolation and in conjunction with all variables). Therefore, for a more conservative interpretation we considered the variables with the top three contributions (i.e., percentage contribution) as important predictors for MaxEnt solution (see Results), and for discussion and interpretation, discarded the rest of the variables due to their low contribution and/or correlation with the three top-performing variables.

We configured the MaxEnt parameters to model the suitable environmental space of each surrogate species (Guillera-Arroita, Lahoz-Monfort, and Elith 2014; Morales, Fernández, and Baca-González 2017). Thus, we obtained logistic values as a proxy of habitat suitability (a map where every pixel has a value between 0 and 1, where values close to 1 represent higher habitat suitability). Furthermore, we used the default value (1) of the regularization multiplier to standardize the error limit of predictive values to avoid their underestimation (<1 more localized distribution) or overestimation (>1 more stretched distribution). Because our occurrence data may have more than one occurrence record for a given pixel, by default the control of duplicate records was activated (see Table 1 for the number of occurrences, training records, and testing records per species). We also selected the default setting for feature types (auto features). Auto features allows the MaxEnt algorithm to automatically select the combination of feature classes used according to the number of species records being modeled as well as standardize the feature selection method (i.e., MaxEnt path) and the response curves of the predictors (environmental data).

We used a 25 percent random test sample of records to evaluate the performance of the MaxEnt models using receiver operating characteristic (ROC) curves, and we explored the area under the curve (AUC) values (Phillips 2017). The main characteristic of ROC analysis is the AUC provides a single measure that is independent of any particular threshold (Phillips, Anderson, and Schapire 2006). We used presence-only data, because the fractional predicted area (the fraction of the total study region predicted present) was used instead of the commission rate (the fraction of absences predicted present). This is particularly important because our approach is directly related to páramo landscapes inside of MCB and our occurrences cover large areas across the entire páramo of the study area, not limited to a few páramo localities (Figure 2). Thus, the maximum achievable AUC value is less than 1 and values close to 1 indicate a good performance of the model because they represent the true positive rate, whereas values around 0.5 are considered no better than a random predictor model (Phillips, Anderson, and Schapire 2006; Morales, Fernández, and Baca-González 2017). Furthermore, AUC maximizes the sensitivity and specificity of ROC, resulting in a good balance between false-positive and false-negative rates (commission/omission) and therefore providing a standard evaluation method (Costa et al. 2015). In addition, we used the omission rate across the proportional predicted area (test point falls into an unpredicted pixel) as an additional method to evaluate how well the model performs. We followed Phillips, Anderson, and Schapire (2006) and we only considered models with less than 5 percent training omission rates.

The models for all eighteen study species showed AUC values >0.96 (mean = 0.99 ± 0.003 SD; Table 1) and all eighteen omission rates were close to zero values (0 percent), indicating good performance. Consequently, we selected the tenth percentile (Table 1) and converted each model into a binary map (predicted/unpredicted). We used the tenth percentile following standard practice (Morales, Fernández, and Baca-González 2017) to produce more interpretable binary maps (i.e., 90 percent of the points were contained within the predicted area). We reclassified the binary maps by assigning the value 1 to grid cells according to the tenth percentile of the output model of each study species (tenth percentile or above = 1, all other values = 0; Figure 2). We then created a stacked map that illustrates the number of study species per pixel (species aggregation) across the study area (Figure 2). The stacked map does not attempt to provide a description of species richness or a gradient of species across the study area (e.g., Calabrese et al. 2014). Rather, this aggregation is a graphical technique (overlaying maps one above the other) to allow the visualization of the intersection of predicted areas (habitat suitability) for all eighteen bird species per predicted pixel. Within this framework, we summed all of the reclassified binary maps together to generate a consolidated map (bring together the models into a single model map) that retained only pixels where all eighteen surrogate bird species are predicted to occur (pixel value = 18; Figure 3). Consequently, the consolidated map is a methodological approach, based on the suitable environmental space, to obtain an overview of where all study species might be present (e.g., key areas for the surrogate bird species) because it is a representation of pixels that contain environmental suitability for all surrogate bird species.

Results

We recorded 2,520 individuals of our eighteen surrogate bird species (mean = 140), with records distributed across both flanks of the study area (Figure 2). Records in the protected areas of MCB represented 43 percent (N = 1,088). The three species with the highest records were many-striped canastero (Asthenes flammulata) with 16 percent of occurrences, tawny antpitta (Grallaria quitensis) with 13 percent, and chestnut-winged cinclodes (Cinclodes albidiventris) with 11 percent of records.

In general, elevation is the environmental variable with the highest contribution to the eighteen MaxEnt models (range = 49–62 percent), followed by precipitation seasonality (range = 15–21 percent) and isothermality (range = 9–15 percent; Table 1). The potential distribution areas (based on the tenth percentile) of the eighteen surrogate species ranged from 115,500 to 312,700 ha (mean = 192,278 ± 43,033.1 SD; Table 1, Figure 2). The three species with the largest predicted distribution areas were shining sunbeam (Aglaeactis cupripennis; total area = 312,700 ha, outside protected area = 88 percent), Ecuadorian hillstar (Oreotrochilus chimborazo; total area = 222,500 ha, outside protected area = 83 percent), and chestnut-winged cinclodes (total area = 217,700 ha, outside protected area = 82 percent). The predicted areas within the system of protected areas of MCB ranged from 30,900 to 38,500 ha (mean = 37,050 ± 184.8 SD; Table 1, Figure 2). The three species with the highest proportion of predicted areas within protected areas were blue-mantled thornbill (Chalcostigma stanleyi; 29 percent), tit-like dacnis (27 percent), and red-rumped bush-tyrant (Cnemarchus erythropygius; 24 percent).

The consolidated map (key areas) covers an area of 94,800 ha and extends from the highlands located southeast of Cañar River in the north, through the mountain systems of Cajas National Park and Quimsacocha National Recreation Area, to the headwaters of Jubones River in the southeast of the study area (mainly between 3,600 and 4,200 m.a.s.l.; Figure 3); 30 percent of this consolidated area (28,200 ha) falls within the protected areas, and the remaining 70 percent (66,660 ha) is outside of the protected areas.

Discussion

The consolidated map of overlapping areas of estimated habitat suitability for eighteen surrogate species suggests the potential conservation importance of the MCB in providing crucial habitat for these bird species. The national system of protected areas provides core of habitat suitability. However, 70 percent of this suitable habitat lies outside of currently protected areas, highlighting additional territories as possible key areas for conservation of high Andean bird species. For example, private lands that are not fully dedicated to livestock and agricultural activities, as well as less accessible public lands, should be considered of conservation value. The conservation of specialized, endemic, and threatened bird species could be improved by the inclusion of these identified areas into more formal protection, such as the categories of the national system of protected areas.

Elevation and precipitation seasonality are environmental variables with important contributions to our MaxEnt models. Our monitoring efforts included both the western and the eastern flanks of the high Andean cordillera; however, the consolidated model predicted greater areas of overlap for surrogate bird species on the eastern flanks, especially at elevational ranges between 3,700 and 4,000 m. This elevational range is crucial because it represents the zone with high availability of páramo habitats for birds across the region (Astudillo et al. 2018, 2019; Barros et al. 2020). This is likely because the prevalent climatic and topographic (especially elevation) conditions on the eastern flank of the study area influence habitat suitability for páramo-specialized birds. In fact, precipitation and elevation are both environmental factors that promote high habitat heterogeneity for páramo habitats, resulting in greater habitat suitability for habitat-specialized organisms across the northern Andes province (e.g., Roach, Urbina-Cardona, and Lacher 2020; Caballero-Villalobos et al. 2021). In addition, at elevations between 3,400 and 3,700 m, the habitats are characterized by a transitional páramo influenced by the tree line of montane cloud forest (Neill 1999). These lower elevations are also subject to pressure from human activities (in particular outside of protected areas) such as grazing (Astudillo et al. 2018) and are particularly sensitive to global change, because some widespread organisms (plants and animals) such as generalist species with rapid dispersal could promote biotic homogenization (McKinney and Lockwood 1999; Foster 2001). As a result, we suggest that a more informative interpretation of habitat suitability across the transitional páramo in particular could be obtained through a widespread monitoring approach in the context of both the entire páramo ecosystem (e.g., provincial páramo) as well as transitional montane páramo habitats, because these páramo habitats in this elevational zone are less available (Astudillo et al. 2018; Barros et al. 2020). In general, the transitional high Andean habitats are influenced by intensification of human activities (e.g., the agricultural frontier) and therefore more precise conservation scenarios are needed (Sarmiento 2000; Sarmiento and Frolich 2002). As seen here, mapping important areas within specific regions may provide an assertive approach to enhancing conservation that is more compatible to a particular ecosystem such as páramo.

On the east-windward slopes (more influenced by the Amazon precipitation regime), there is relatively higher rainfall with less seasonality that leads to greater cover of native tussock grasslands associated with cushion plants and woody plants in humid areas (Sklenár and Ramsay 2001; Ochoa-Sánchez, Crespo, and Célleri 2018; Lazo et al. 2019). The influence of high rainfall with less seasonal regimes leads to more habitat availability for grassland-specialized birds (e.g., Antos and Schultz 2020). Our individual models show that species that use or are specialized in páramo grassland habitats, such as violet-throated metaltail, Andean tit-spinetail (Leptasthenura andicola), stout-billed cinclodes, and tit-like dacnis, are important components of the bird community. A few studies within the protected areas of the MCB previously reported that the diversity of páramo-specialized birds increased with páramo heterogeneity (Tinoco et al. 2009; Astudillo et al. 2019, 2020), and this is commonly reported at relatively higher elevations (Astudillo 2014). Here, our study not only corroborates these observations but maps out important areas where suitable habitat for páramo-specialized birds may be found, especially outside of the protected areas. In particular, unprotected páramos are located beyond the northeastern limit of PNC, the eastern flank between PNC and ARQ, as well as beyond the southeastern limit of ARQ.

On the western flank (facing the Pacific basin), precipitation seasonality shows a more bimodal regime with two events, especially compared with similar elevational ranges on the eastern flank (Campozano et al. 2016, 2018). This particular climatic regime leads to less tussock grass habitat, with drier areas and increasing cover of woody plants (Neill 1999; Koenen 2000; Barros et al. 2020). This results in lower habitat suitability for more grass-dependent páramo birds. Upon investigating individual species responses, we found that some hummingbird species such as Ecuadorian hillstar and shining sunbeam and songbirds such as many-striped canastero and brown-backed chat-tyrant (Ochthoeca fumicolor) have an important predicted suitable area on this western flank. These species have greater affinity to semi-open habitats with higher cover of woody bushes and shrubs (Astudillo, Tinoco, and Siddons 2015; Barros et al. 2020). These birds are also key species for providing ecosystem services such as pollination and insect control (Astudillo et al. 2019; Cárdenas et al. 2020). Our study suggests a number of areas where habitat characteristics indicate importance for these bird species, including surrounding páramos located beyond the northwestern and western limits of PNC and beyond the southwestern limit of ARQ. Managers of natural areas may be guided toward conservation decisions by recognizing more localized areas as important elements to be considered.

Conservation remarks

The MCB and its protected areas have been recognized as a priority region for biodiversity conservation in Ecuador (e.g., Lessmann, Muñoz, and Bonaccorso 2014; Cuesta et al. 2017). However, this recognition is broad and not specific for the páramo grassland landscape. For instance, proposed extensions to Cajas National Park by Lessmann, Muñoz, and Bonaccorso (2014) include areas to the north identified in this study and inclusion of montane forest on the western slopes. These identified areas do not include páramos in southern MCB or around ARQ as a priority; consequently, these extensions do not consider the highest mountain habitats. Highest habitats in the Andean region are a priority for conservation because they relate to remarkable biodiversity associated with high-quality ecosystem services (e.g., water regulation, carbon storage; Sarmiento 2000; Carrillo-Rojas et al. 2019; Lazo et al. 2019). Our findings highlight páramo areas worthy of consideration for avian protection outside of present protected areas because 70 percent of the consolidated model falls outside of these areas. In particular, the eastern flank of the MCB shows suitable habitat for grass-specialized bird species, especially for birds that prefer more open tussock grass habitats with less woody native plants. In addition, the western flank of the MCB should be considered a priority for localized conservation plans for birds specialized in shrubby páramo habitat. This study improves our understanding of vital areas of páramo birds via habitat suitability, especially in relation to the environmental predictors, and contributes critical data for setting conservation priorities in high-altitude habitats associated with these bird species within the páramo ecosystem (the highest vegetated frontier in the northern Andes). These data can be used to identify new potential areas for protection within the regional páramo. Though this research is limited to the study region due to the uniqueness of the biogeography of the MCB (Jiménez-Rivillas et al. 2018), the approach demonstrated here, based on systematic surveys across the high Andes and species modeling, is applicable regionally for northern Andes grasslands with high climatic variation (Lees et al. 2020; Soares et al. 2023).

Furthermore, findings from this study highlight the general importance of areas currently outside of the national system of protected areas for consolidating high Andean bird conservation. This knowledge is important because habitat loss through land use change (i.e., livestock grazing, road infrastructure, exotic tree plantations) is widespread across the study region (Astudillo 2014; Aguilar et al. 2019). On a landscape scale, there is a large body of evidence that these human activities are associated with decreasing availability of páramo habitats (e.g., Sarmiento 2000; Suarez and Medina 2001; Latta et al. 2011; Matson and Bart 2013). This results in negative impacts on ecosystem functionality (Matson and Bart 2014), leading to decreasing diversity of both grass-dependent páramo-specialized birds and shrubby paramo–specialized birds (Astudillo et al. 2018; Aguilar et al. 2019). Suitable habitat for páramo birds is not only found in protected areas. For instance, the approach of spatial differentiation of birds (e.g., regionalization of biodiversity concentration) in the northern Andes is postulated as an important factor to locate potential territories to designate as protected areas (e.g., via habitats that contain populations of special-interest species or surrogate species; e.g., Kattan et al. 2006). This study shows that there are potential priority areas for conservation located on private lands, especially those that are not fully dedicated to livestock or agricultural activities (Barros et al. 2020), as well as public lands that are less accessible to human activities at relatively higher elevations. On such lands, the habitat heterogeneity of páramo (i.e., tussock grass in combination with native woody plants) harbors a comparatively high diversity of paramo-specialized birds (Astudillo et al. 2018, 2019; Barros et al. 2020). Consequently, findings from this study represent an effort, based on ecological modeling, to improve conservation policies across the region. By overlaying the consolidated map over polygons of public and private lands with little evidence of changes in land use, areas with higher habitat suitability for bird conservation can be identified and therefore guide conservation toward Andean regions with unique diversity and high endemism, such as in the Macizo del Cajas Biosphere Reserve.

Authors’ contribution

PXA: Conceptualization, methodology, research, resources, field data, writing—preparation of the original draft, supervision, acquisition of funds, project fund acquisition, project management. SB: Methodology, field data, writing—preparation of the original draft. DM: Conceptualization, methodology, writing—preparation of the original draft. FRV: Methodology, writing—preparation of the original draft. DCS: Methodology, writing—preparation of the original draft. SCL: Conceptualization, methodology, writing—preparation of the original draft, supervision, acquisition of funds.

Ethics/legal

The authors declare that they did not violate or omit ethical or legal norms in this research. All activities were conducted under permit numbers 196-2019-DPAA/MA and MAAE-ARSFC-2021-1075 issued by the Environmental and Water Ministry of Ecuador. The study data are freely available from the corresponding author.

Acknowledgments

We are grateful for the consistent support of our research activities by Jacinto Guillén, Raffaella Ansaloni, and Andrés López from the Universidad del Azuay. We thank Boris Landázuri, Eduardo Barnuevo, Paúl Porras, Bruno Timbe, and Ramiro Jiménez for their field assistance. We also acknowledge Nikolay Hristov, Scott Campbell, Rafael Salgado, and Vicente Jaramillo from Dundee Precious Metals Ecuador (DPMECUADOR S.A.) for their continued support of our research. Logistical support came from Hari Gonzáles, park management, and staff of Cajas National Park and Quimsacocha National Recreation Area, as well as Empresa Pública de Telecomunicaciones, Agua Potable, and Alcantarillado y Saneamiento de Cuenca (ETAPA-EP). This article greatly benefited from comments and editing by Omar Delgado, Gary Ritchison, Reed Bowman, and Nina Attias.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15230430.2023.2299362

Additional information

Funding

PXA was funded by Vicerrectorado de Investigaciones, Universidad del Azuay [Grant Nos. 2017-139; 2018-0003; 2020-0174] and ETAPA-EP [Grant No. 78010200001]. SCL was funded by an Avian Conservation Endowment of the National Aviary (USA). The field activities were funded by INVMETALS (now DPMECUADOR S.A.).

Appendix

Correlations of six environmental variables used to predict, via MaxEnt algorithm, the habitat suitability of eighteen surrogate species in the páramo landscape of Macizo del Cajas Biosphere Reserve, southern Andes of Ecuador. The codes for bioclimatic variables (BIO) are as follows: BIO18, precipitation amount of warmest quarter (kg mm⁻²); BIO15, coefficient of variation of precipitation (precipitation seasonality); BIO16, precipitation amount of wettest quarter (kg m⁻²); BIO2, mean temperature of diurnal range (°C); BIO8, mean temperature of wettest quarter (°C); and BIO3, isothermality (relation between the diurnal temperature range and annual temperature range). These environmental variables were selected because they represent a relatively high contribution to previous MaxEnt outputs at the scale of páramo-specialized species across the study area (Tinoco et al. 2009; Astudillo et al. 2016). The influence of correlated variables to the MaxEnt optimal solution was evaluated by the importance of percentage contribution of models as well as by jackknife checkbox (see Methods).