Bird species richness in High-Andean forest fragments: habitat quality and topography matter

1 Evolutionary Ecology Group, Department of Biology, University of Antwerp, Universiteitsplein 1, 2610 Wilrijk, Belgium. 2 Terrestrial Ecology Unit (TEREC), Ghent University, K.L. Ledeganckstraat 35, 9000 Gent, Belgium. 3 Center for Macroecology Evolution and Climate, Universitetsparken 15, University of Copenhagen, 2100 Copenhagen, Denmark. 4 Centro de Biodiversidad y Genética, Facultad de Ciencias y Tecnología, Universidad Mayor de San Simón, Calle Sucre frente al Parque La Torre, Cochabamba, Bolivia.


Introduction
Montane forests worldwide are important centers of biodiversity and endemism and provide numerous ecosystem services as they help to retain water, stabilize soils, store carbon and increase soil fertility (aubaD et al. 2010). Unfortunately, montane forests, located in topographically complex areas, are heavily affected and threatened by fragmentation processes, mainly caused by anthropogenic activities such as land conversion for agriculture (Doumenge et al. 1995;Cayuela et al. 2006). To protect and restore these habitats and their biodiversity, it is important to understand factors and processes influencing species persistence (Flaspohler et al. 2010;tinoCo et al. 2013). Such information is especially urgent for the management of Polylepis woodlands, a unique ecosystem with large numbers of endemic and threatened species occurring in the South American High Andes up to 5200 m (FJelDså 1993;navarro et al. 2005). Due to extensive historical and ongoing anthropogenic activities, less than 10% of the original Polylepis forest cover is estimated to remain, making it one of the most endangered forested ecosystems in the world (WCMC 2004). This degradation is mainly caused by slash-and-burn agriculture techniques, cattle overgrazing, firewood collection and the replacement of native forests by exotic plantations (FJelDså & Kessler 1996;hensen 2002;balDerrama 2006;gareCa et al. 2007gareCa et al. , 2010hensen et al. 2012). Today, Polylepis forests mostly remain as small patches restricted to highly inaccessible areas like ravines, ledges and steep slopes (renison et al. 2011; sylvester et al. 2014; alinari et al. 2015). There is a strong debate as to what extent their patchy distribution is natural (gosling et al. 2009;Williams et al. 2011) or caused by anthropogenic activities, but evidence points to human-induced pressures leading to extensive habitat loss and increased isolation over time (FJelDså & Kessler 1996;CierJaCKs et al. 2007;torres et al. 2008;toivonen et al. 2011;alinari et al. 2015).
While little is known of the characteristics of Polylepis forest patches influencing biodiversity within Bolivia, the area, elevation and habitat quality of forest fragments all have been found to influence Polylepis bird communities in forests of Argentina, Peru and Ecuador (FJelDså 1993;lloyD 2008a;bellis et al. 2009, 2015sevillano ríos et al. 2011;tinoCo et al. 2013;sevillano-ríos & roDeWalD 2017). Although these studies showed that small forest fragments may constitute important habitats for many species, including threatened and endemic ones (lloyD & marsDen 2011), larger and/or higher quality forest remnants were often found to support larger and more diverse communities (FJelDså 1993;bellis et al. 2009, 2015. Predictors such as soil erosion, humidity, habitat complexity, Polylepis tree species and plant composition can also affect Polylepis bird communities (FJelDså 1993;lloyD & marsDen 2008;bellis et al. 2009, 2015tinoCo et al. 2013;sevillano-ríos & roDeWalD 2017). However, few fine-scale topographical features, beside elevation and slope angle, have been included in such analyses ( lloyD & marsDen 2008;sevillano ríos et al. 2011;bellis et al. 2015) despite their role in increasing habitat diversity locally (homeier et al. 2010).
Here, we assessed which factors correlate with bird species richness and community composition in a fragmented Polylepis forest area in the Tunari National Park (TNP), one of the most important areas for the conservation of Polylepis and its associated avifauna in Bolivia and South America in general (FJelDså 2002;balDerrama 2006;Faunagua 2015). The park was created in 1962 to protect the native vegetation, to improve water resource management and to prevent erosion, landslides and floods. It should also limit the expansion of the city of Cochabamba and function as a recreational area for its residents, yet urban encroachment and land conversion to agriculture have taken a heavy toll, particularly on its Cochabamba-facing southern slope (sernap 2016). Nevertheless, the TNP still represents the last main stronghold of two vulnerable tree species endemic to Bolivia, namely Polylepis subtusalbida (distributed on the southern slope) and P. lanata (mainly present on the more humid northern and eastern slopes). Both forest types are important to support endemic and threatened animal and plant species (FJelDså 2002;gareCa et al. 2010). The southern slope area is specifically important for the conservation of the endangered Cochabamba-mountain finch (Compsospiza garleppi, IUCN 2018) and several High-Andean species and has consequently been designated an Important Bird and Biodiversity Area by BirdLife International (2017).
Because in the TNP, Polylepis patches are strongly clustered in different watersheds, we investigated factors driving bird species richness patterns and community composition across different watersheds ('watershed-level analysis') before focusing on individual patches ('patch-level analysis'). In addition, we studied which factors affect the presence of endemic and/or threatened bird species at both watershed and patch levels. We hypothesize that factors related to fragmented or poor habitat quality patches, such as smaller fragment size, isolated and small Polylepis trees and the presence of exotic plantations, will negatively affect bird species richness and bird community composition, in particular the presence of threatened and specialized species. Inversely, topographical complexity of fragments with large smallscale terrain heterogeneity and (topography-driven) soil moisture is hypothesized to allow for more diverse bird communities.

Study area
The Tunari National Park (TNP, 16°55′-17°34′ S, 66°55′-66°44′ W), a 300 000 ha area located in Bolivia, ranges from 2750 to 4400 meters a.s.l. Above 3200 meters, two endemic vegetation associations dominate in the park (navarro et al. 2005): Berberis commutata -Polylepis subtusalbida and Styloceras columnare -P. lanata (Kessler & sChmiDt-lebuhn 2006). The first woodland type is characteristic of the arid habitats found on the southern slope of the park (precipitation: 600-800 mm) and forms monospecific tree stands while the second is associated with trees of the yungas habitats on the more humid northern and eastern slopes (precipitation: 900-1000 mm, navarro et al. 2005). Native forests are estimated to cover about 4% of the park, although the exact distribution of the remaining Polylepis patches remains unknown (sanabria siles et al. 2012).

Forest mapping
Polylepis remnants were located in the TNP with a combination of available maps and a LANDSAT8 image (30 m resolution). Their exact boundaries were mapped with a GPS device during field surveys conducted between September and December in 2014 and 2015 and complemented with Google Earth for inaccessible areas. These remnants are generally strongly isolated as they are restricted to watersheds located in different valleys separated by deep ravines or high ridges. On the southern slope, remnants within watersheds are further fragmented into smaller patches which are separated from each other by agricultural fields, pastures, puna grasslands or exotic plantations (Eucalyptus globulus and/or Pinus radiata). Because of this clustered configuration of the patches, we decided to carry out the analysis at 'watershed' and at 'patch' levels, assigning each Polylepis patch to a specific watershed. Some watersheds only contain a single patch (Fig. 1, Table 1). We uncovered four watersheds containing Polylepis lanata patches (one patch per watershed, totaling some 1252 ha) and ten watersheds with Polylepis subtusalbida (one to three patches per watershed, in total 16 different patches, which together cover about 832 ha; details in Tables 1-2). Our patch-level analysis considers all sampled individual patches of Polylepis subtusalbida of the southern slope.

Bird surveys
Bird surveys were carried out from September to December in 2014 and 2015. We surveyed nine out of ten watersheds where P. subtusalbida patches were located on the southern slope (corresponding to 15 of the 16 P. subtusalbida patches uncovered) and two P. lanata watersheds (Fig. 1). We were unable to survey the remaining three watersheds (one P. subtusalbida and two P. lanata watersheds), for which we relied on previously collected data (see below). We used an adapted version of the MacKinnon-list method, which groups observations of birds into consecutive lists of species. A species accumulation curve is generated by adding those species not recorded in any previous list to the total species number, which is then plotted as a function of the list number. This method thus relates cumulative species Fig. 1 -Map of the Tunari National Park (Bolivia), showing elevation (in meters) and all Polylepis forest fragments included in the analysis. Names of the watersheds are indicated in bold outside the frames, in black for Polylepis subtusalbida and in grey for Polylepis lanata tree species. Patches within watersheds are shown within frames. * indicates watersheds for which we obtained bird data from a local expert.  1 Polylepis forest watersheds inside the Tunari National Park with the number of patches and total size of forest cover detected and surveyed per cluster, Polylepis tree species, observed species richness (Sobs), estimated species richness (Chao2 and ICE), the number of species observed by José Balderrama ("Expert list", used when no other data was available) and the total survey time in minutes.
richness to the number of observations, rather than time or space, and thereby accounts for moderate differences in observer qualification and field conditions (maCKinnon & phillips 1993(maCKinnon & phillips ). herzog et al. (2002 found that for tropical, species-poor habitats (such as Polylepis remnants), 10-species lists were more robust. Therefore, our master list of observations was later processed into such 10-species lists, each list starting with the observation following the last observation of the previous list and therefore containing the ten first different species observed (herzog et al. 2002;o'Dea et al. 2004;maCleoD et al. 2011;Cavarzere et al. 2012). In addition, following herzog et al. (2002), observations beyond 50 m from the observer were excluded as detectability of forest birds substantially declines beyond that distance.
Surveys were carried out by two skilled observers (C.F and H.L, the latter replaced by J.B for the surveys conducted in MOR, Fig. 1). Both observers walked slowly and randomly inside the patches and recorded the species (identified visually or orally) for each encounter and the time of observation. Unknown songs and calls estimated to originate within a 50 meter-radius from the observer were recorded to be later identified and incorporated into the species lists. As they are difficult to differentiate in the field, Turdus chiguanco and T. fuscater were both recorded as Turdus spp. (TURSP) and Nothoprocta pentlandii and N. ornata were both recorded as Nothoprocta spp. (NOTSP). Surveys were conducted on days without wind or rain from dawn (05h30) to midday and from late afternoon to dusk (18h30) for a total of 321.6 hours. Polylepis patch size ranged from 3.2 to 1028.9 ha (median: 17.9 ha, Table 2 and Figs S1-S2), and survey intensity was roughly proportional to patch size (r = 0.50, P = 0.055, Table 2). Flyovers and obvious non-target species (not specifically associated with forest habitats, e.g., water-birds, grassland specialists and raptors) were excluded. For the three watersheds that we could not survey ourselves (i.e., LIR, ICA and ESP; Fig. 1  these patches to sample all bird species present, during which J.B. walked slowly across the forest patches, actively searching any bird species previously undetected. Thus, our final dataset consists of bird surveys in ten watersheds located on the southern slope (totaling 16 P. subtusalbida patches), and four watersheds elsewhere in the park (each of them containing one P. lanata patch).

Vegetation characteristics
For the 15 P. subtusalbida patches where we conducted bird surveys ourselves, vegetation characteristics were obtained in 279 randomly selected 5 m diameter-sampling plots (approximate density of one plot per ha). In these, we visually estimated (1) the proportion of exotic trees within the overall tree cover ('exotic tree ratio') and (2) tree height, and we measured (3) diameter at breast height (DBH) and (4) canopy density with a convex spherical crown densimeter. We scored (5) the degradation level caused by fire, logging and grazing on a scale from 0 to 3 (none, low, medium or high) and used the median of these scores as an 'anthropogenic degradation index'. Due to difficulty of access, these measures are not available in three of the watersheds (LIR, ICA and ESP, Fig. 1) and were thus not included in the watershed analysis. Lastly, we mapped all exotic plantations (Eucalyptus globulus and/or Pinus radiata stands) located within 150 m of the Polylepis remnants with a GPS device. Because between-patch movements of Polylepis-associated birds decrease at distances of more than 200 meters (lloyD & marsDen 2011), we calculated (6) the extent to which each Polylepis patch/watershed is surrounded by these plantations as the average proportion of the surface covered by exotics across three buffer areas at 50, 100 and 150 m distances (adapted method from DunForD & FreemarK 2005).

Spatial habitat variables
Using a Digital Elevation Model (ASTER GDEM, 30 m resolution, https://lpdaac.usgs.gov/) in ArcGIS 10.1, we derived for each Polylepis patch (7) mean elevation, (8) (log-transformed) surface area, (9) the average number of hours of sunlight at least half a watershed or patch receives per day ('illumination time'), (10-13) four different measures of local topographic variation (percentage of depressions, ridges, steep slopes and gentle slopes of each watershed or patch) derived from the Topographic Index (TPI, see S1 for detailed technical explanation), (14) the Topographic Wetness Index, (TWI, beven & KirKby 1979) a measure of the topographical control over hydrological processes which is related to forest productivity and soil moisture (besnarD et al. 2013;Wilson et al. 2013;Campos et al. 2015). The TWI represents the potential accumulation capacity in a given pixel and was obtained through the D-infinity flow accumulation algorithm with the TauDEM tool (tarboton 2009) in ArcGIS 10.1. Finally, we computed from LANDSAT8 images (15) the mean Normalized Difference Vegetation Index (NDVI), a proxy for resources available to consumers (hurlbert & hasKell 2003;bellis et al. 2015) and (16) the Tasseled Cap Wetness index (TCW), the actual soil and plant moisture content, with the Tasseled Cap function in ArcGIS 10.1. LANDSAT images (30 meters resolution, courtesy of the U.S. Geological Survey) were recorded in August 2015, the driest month of the year, when the contrast between dry puna grasses/agricultural crops and evergreen Polylepis (and other native) trees or exotic plantations is maximal. All calculations were performed at watershed and patch levels separately.
To avoid the inclusion of correlated variables, two Principal Component Analyses (PCA) were performed for the patch-level analysis. A first PCA ('TopoPCA') included all variables describing patch topography (variables (7), (9), (10-13)) while a second PCA ('VegPCA') was performed on habitat quality and vegetation-related variables ((1), (2), (3), (4), (5)) as well as (6) surrounding exotics and (15) mean NDVI. At watershed level, one PCA was applied for topographic variables while surrounding exotics and NDVI were included as separate variables due to the lack of vegetation data. PCA was carried out with the package 'ade4' in R (Dray & DuFour 2007). We retained the number of PCA axes needed to capture > 80% of the variation present in the data as explanatory variables (three axes at the patch-level, two at the watershed level, Tables 3-4). As they are uncorrelated to other variables, area and both moisture-related variables (TWI and TCW) were used as separate explanatory variables. At both patch-and watershed-level, no strong correlations remained between final sets of explanatory variables (patch-level: all r < 0.48, watershed level: all r < 0.67). At patch level, the first axis of the TopoPCA (TopoAxis1) represents well illuminated patches with a smooth topography. The second axis (TopoAxis2) is associated with high elevation and rugged (i.e., high numbers of depressions in the terrain) patches, while the third axis (TopoAxis3) is correlated with high elevation, high illumination, and rugged terrain (prevalence of steep slopes and high numbers of depressions). At watershed level, we retained as explanatory variables the two first axes of the TopoPCA (TopoAxis1 and TopoAxis2), with similar interpretations as above. The first axis of the VegPCA (VegPCA1) represents forest patches strongly affected by exotic plantations, which contain and are surrounded by high amounts of exotics with sparse, small Polylepis trees. The second axis (VegAxis2) is indicative of high-quality patches, i.e., denser patches with tall Polylepis trees and limited anthropogenic degradation levels. The third axis (VegAxis3) is correlated with patches heavily affected by exotic plantations, characterized by little and small Polylepis trees, with many exotic plantations inside or around the patches.

Bird species richness estimates
Because our species accumulation curves show that observed species richness does not reach an asymptote in the sampled watersheds and patches (

Statistical analysis
To explain species richness patterns at watershed and patch level, we applied Generalized Linear (Mixed) Models (GL(M)M) in which Chao2 or ICE species richness estimates were specified as dependent variable. At the watershed-level, fixed effects included the first two axes of the watershedlevel topographical PCA, (log) area, NDVI, surrounding exotics, TWI, TCW and finally Polylepis tree species. At the patch-level, in addition to the first three axes of the topographical and vegetation PCAs, (log) area, TWI and TCW were modelled as fixed effects while the watershed in which a patch is situated was included as a random effect. Following burnham & anDerson (2002), model selection proceeded through a multi-model averaging framework based on the Akaike's Information Criterion corrected for small sample sizes (AIC C ), as implemented in the R package MuMin (barton et al 2018). We ran all possible combinations of predictor variables, and as no single model attained decisive support (i.e., AIC C weight (AIC W ) > 0.95), AIC W were calculated for all variables (Table S1). Only variables with AIC W > 0.5 were considered informative (burnham & anDerson 2002). We validated the models by using random subsets of the watersheds/patches with a bootstrap analysis of 100 random partitions in R (r Development Core team 2008). Normality of model residuals was verified and met (Shapiro-Wilk W > 0.9). All statistical analyses were performed with the software R.
A separate analysis was conducted for each of the threatened and/or endemic bird species recorded more than once in the area (based on IUCN 2018, Table S2), at both watershed and patch level. These were five threatened species (Asthenes heterura (near threatened, NT), Compsospiza garleppi (endangered, EN), Sylviorthorhynchus yanacensis (NT), Conirostrum binghami (NT), Pseudosaltator rufiventris (NT) and four country endemics (Oreopsar bolivianus, Aglaeactis pamela, Coeligena violifer, Asthenes harterti). We applied a similar modelling framework as described above, but because we focused on the presence/absence of these species in different watersheds or patches, a binomial error distribution was used. To avoid model fitting difficulties related to (quasi) separation, we applied Fisher's penalized logistic regression as implemented in the R package logistf (heinze et al. 2013). As above, predictors were considered as potentially relevant when their AIC W were > 0.5.
To investigate community composition patterns at watershed and patch level, we applied non-metric multidimensional scaling (NMDS; bray distance) based on bird species presence/absence data in watersheds and patches, respectively (Table S3). Generally, ordination with stress values smaller than 0.2, 0.1 and 0.05 are considered of fair, good and excellent quality, respectively. We then identified which species most strongly underlie the ordination (p < 0.05) and which environmental factors (the same factors used to explain species richness, see above) correlate with the ordination axes (p < 0.05) for each axis of the NMDS with the envfit function from the R package vegan (oKsanen et al. 2019).

Bird species richness
We recorded 144 species, of which 50 were found exclusively in Polylepis lanata patches, and 80 only in P. subtusalbida (Table S3), among which the endangered and endemic Compsospiza garleppi. At the watershed level, Chao2 estimates of species richness vary from 33.5 (CHO) to 71.7 (TAQ), with a mean of 48.5. ICE estimates vary from 34.6 (CHO) to 70.2 (TAQ), with a mean of 47 (Table 1). At the patch level, observed species richness varies from 22 bird species (SAM3) to 44 (TAQ). Chao2 estimates ranged from 30.9 (PINT1) to 71.7 (TAQ), with an average of 48.2. ICE estimates ranged from 32.0 (PINT1) to 70.2 (TAQ) with an average of 42.8 (Table 2).

Species richness predictors
Different variables influenced overall bird species richness at the watershed versus patch level. Topography was most important at the watershed level, as watersheds with a high topographic potential to retain rainwater harbored more species (Chao2 and ICE; AIC W = 0.62 and 0.52, Fig. 2a, Table S1).
At the patch level, species richness was highest in productive Polylepis patches with few exotic plantations, regardless whether Polylepis trees were tall and dense (Chao2 and ICE estimators with AIC W = 0.60 and 0.63, respectively) or not (Chao2 estimator with AIC W = 0.68; Fig. 2b, Table S1). Species richness was also higher in rugged low-elevation patches with short illumination times and few steep slopes (Chao2 and ICE estimators with AIC W = 0.88 and 0.87, respectively). Surprisingly, fragment area did not appear to be an important predictor of bird species richness at either patch or watershed level (i.e., all AIC W < 0.28; Tables S1-S2).

Birds of conservation concern
Four out of the nine bird species of conservation concern were only included in the analysis at the watershed level, either because they were absent from P. subtusalbida patches (Aglaeactis pamela, Coeligena violifer, Asthenes harterti) or present in all of them (Pseudosaltator rufiventris). Topographyrelated variables (TopoPCA axes) influenced the presence of all species of conservation concern (Table 5), except C. violifer, whose distribution was mainly governed by vegetation productivity (high NDVI values). Vegetation parameters, particularly the two first axes of the VegPCA, influenced the patch-level distribution of four species (Table 5). Compsospiza garleppi was only observed in P. subtusalbida watersheds located at higher elevation and on steep slopes, where it was most likely to be observed in patches with small Polylepis trees but unaffected by exotic plantations. The Polylepis specialist Conirostrum binghami occurred more often in the most elevated and rugged P. subtusalbida watersheds, and in the higher elevated, well illuminated and less degraded patches. The other specialist, Sylviorthorhynchus yanacensis was mainly present in well-lit P. subtusalbida stands, particularly in drier and less degraded patches. Asthenes heterura was most likely to be found on well-illuminated, gentle slopes while Oreopsar bolivianus was more observed in high-elevation and steep P. lanata watersheds, devoid of exotic plantations. The occurrence of P. rufiventris was limited to elevated P. subtusalbida watersheds on rugged terrain. A. harterti occurred exclusively in P. lanata stands, avoiding rugged areas at higher elevations. Finally, A. pamela was observed in P. lanata forest watersheds located on wellilluminated gentle slopes.

Community composition
Thirty-three species out of 144 contributed significantly to the NMDS ordination of the watersheds (Table S4) and two factors, Polylepis species and NDVI, were significantly associated with the first axis of the NMDS (p = 0.05 and p = 0.047 respectively, Fig. 3a). None of the tested factors were significantly associated with the second axis of the NMDS. The ordination was good as residual stress value equaled 0.099 (non-metric fit R² = 0.979). Nine bird species, among which the endemic A. harterti, A. pamela and C. violifer, significantly drove the ordination on the left of the first axis, which corresponds to P. lanata watersheds characterized by high NDVI values (Fig. 3a, Table S4). Eight species, among which the near-threatened S. yanacensis and P. rufiventris, significantly drove the ordination to the right of the first axis, which corresponds to less productive P. subtusalbida watersheds (Fig. 3a, Table S4). Fourteen species significantly drove the ordination on the second axis, a pattern which was not related to any of the factors included in the analysis, but upon visual inspection appears to mainly correspond to the P. subtusalbida LIR watershed (Fig. 3a, Table S4).
Twenty-two bird species significantly drove the NMDS ordination of patches on at least one of the NMDS axes (Table S5). Two factors related to habitat quality and topography were significantly associated with the first axis of the NMDS (VegAxis2, p = 0.048 and TopoAxis2, p = 0.02, respectively; Fig. 3b) while one factor related to topography was significantly associated with the second axis (TopoAxis1, p = 0.002). The ordination was fair as residual stress value equaled 0.144 (non-metric fit R² = 0.99). Three species, among which the near-threatened C. binghami and S. yanacensis, significantly drove the ordination on the right of the first axis, which corresponds to denser, less degraded and more elevated patches with rugged terrains while nine species, none of them of conservation concern, significantly drove the ordination on the left side of the first axis, which corresponds to low elevation degraded patches (Fig. 3b, Table S5). Six species, among which the endangered C. garleppi, significantly drove the ordination to the bottom of the second axis, which corresponds to well-illuminated patches located on smoother terrains, while four species without conservation concern significantly drove the ordination to the upper section of the second axis, which corresponds to less illuminated patches located on topographically complex terrain (Fig. 3b, Table S5).

Discussion
Polylepis forest remnants of the Tunari National Park constitute an important habitat for Andean bird conservation as they harbor 144 different species, including numerous endemic and/or threatened species (Table S2). We observed 18 out of the 26 bird species that were originally used to support the designation of the park as an Important Bird and Biodiversity Area, highlighting the important role of small Polylepis remnants for biodiversity. Overall avian species richness, the distribution of species of conservation concern and community composition across the TNP were influenced by both topography and Polylepis remnant habitat quality. Moister Polylepis remnants harbored richer bird communities, while on the southern slope, bird species richness was highest in P. subtusalbida patches located at lower elevations and on less steep but uneven terrains with high amounts of local depressions. The presence of exotic plantations, in-or outside the patches, negatively affected bird species richness, and rugged Polylepis patches characterized by high exposure to sunlight, denser canopies and high cover of tall and thick Polylepis trees were more likely to support species of conservation concern. Polylepis lanata and P. subtusalbida forest remnants generally tended to support different species, including several species of conservation concern.
The topographical position of Polylepis patches influences its associated avifauna through both their local water retention capacity and the amount of sunlight they receive. Across the park, we found that areas located in basins or valleys with a topology favoring rainwater retention, and therefore potentially being characterized by higher soil moisture, exhibit more diverse communities. Water availability is a major limiting factor for tree growth in the high Andes, even though Polylepis trees are physiologically adapted to survive in arid environments (azoCar et al. 2007 On the southern slope, we found that bird communities of well-illuminated P. subtudalbida patches were less diverse but contained distinct communities characterized by the presence of several species, among which the endangered and endemic C. garleppi. A positive association with exposure to sunlight was also found for several species of conservation concern, particularly insectivorous A. heterura, S. yanacensis, C. binghami and nectivorous A. pamela. The amount of sunlight reaching a forest patch characterizes its microclimate and therefore influences its insect activity and plant productivity. Though such a sun exposure effect has already been reported for C. binghami (De Coster et al. 2009), our results suggest that illumination is also important for other bird species in montane environments and influences the composition of communities inhabiting forest patches.
We found more species-rich, distinct bird communities in rugged forest remnants, probably because topographically diverse areas offer more suitable breeding, roosting and sheltering sites, or protection from predators (martínez-morales 2005). Uneven grounds also protect the soil against erosion caused by grazing livestock (torres et al. 2008; renison et al. 2010; alinari et al. 2015; bellis et al. 2015). Topography appears to be especially important for C. binghami and P. rufiventris which are largely restricted to forest stands located at high elevations, on steep slopes and rugged terrains.
Besides topography, we found that forest remnants surrounded by pine or eucalyptus plantations supported less diverse communities, a pattern which has been observed in other . Exactly why the presence of exotic plantations in and around Polylepis fragments negatively affects bird species richness and the presence of some species of conservation concern inside the fragments remains unclear. A first possibility is that they reduce habitat quality within forest remnants. Exotics, and eucalyptus trees especially, can indeed affect soil quality, compete with the native vegetation and become invasive (gareCa et al. 2007;branDt et al. 2012;brugger et al. 2019). A second possibility is that these dense plantations of high trees surrounding forest remnants constitute a barrier to many bird species unable or unwilling to cross large distances in unsuitable habitats (lloyD & marsDen 2011) and therefore increase the level of isolation of already heavily fragmented and scattered Polylepis fragments. Whichever processes are involved, and given our findings, we would recommend to prevent the establishment of exotic plantations in the vicinity of Polylepis fragments in the TNP (mortelliti & linDenmayer 2015).
It should, however, be noted that responses to Polylepis fragmentation and matrix composition likely are species-specific (preveDello & vieira 2010). This is well illustrated by the endangered and endemic Cochabamba mountain-finch (Compsospiza garleppi), whose range is restricted to the mountain slopes in the departments of Cochabamba and Potosí (balDerrama 2009; birDliFe international 2012). Even though long considered a strict Polylepis-specialist, it was recently found to be rather dependent on Polylepis-associated shrubs and was also observed feeding on nearby crops (huanCa et al. 2009). Our results show that Cochabamba mountain-finches are indeed more tolerant to anthropogenic exploitation and forest degradation than true specialists such as C. binghami and S. yanacensis (although C. garleppi is less likely to occur near forest fragments containing or being surrounded by exotic plantations). In the TNP, C. garleppi however has a privileged relationship with the P. subtusalbida ecosystem as it was only observed in such fragments, particularly in well illuminated and rugged patches that were not surrounded by exotic plantations.
Finally, but importantly, our results indicate that patch surface area did not serve as an important predictor for bird species richness, nor for the presence of any of the species of conservation concern. Reported effects of fragment area on species richness inside Polylepis patches have been mixed so far, with some studies reporting a strong negative correlation (tinoCo et al. 2013) while others found none (lloyD 2008a). FJelDså & Kessler (1996) also observed that the presence of Polylepis-associated bird species across the Andes was not strongly related to patch area and that many species persisted in most of the remaining patches, if they were of sufficient quality. It is thus generally assumed that Polylepis-dependent species have become adapted to the fragmentation of their habitats due to the long history of anthropogenic deforestation in the area and that they survive as metapopulations, i.e., small populations connected by dispersal, although some species have been reported to avoid tiny (< 1 ha) isolated patches (lloyD & marsDen 2011). Unfortunately, little is known about the number, size and spatial configuration of the fragments necessary for the metapopulations of these species to persist (hansKi 1998;hansKi & ovasKainen 2002). More research is therefore needed to know if the small populations recorded in the small Polylepis patches of the TNP are viable or if they are expected to go extinct (extinction debt, purCell et al. 2002;hansKi & ovasKainen 2002;lloyD 2008b;sevillano ríos et al. 2011). Understanding such processes is particularly important for specialized and threatened species which occur at very low densities in the area such as C. binghami and L. yanacensis (Cahill & matthysen 2007).
There were several logistical constraints associated with surveying the avifauna in the remote and heavily fragmented montane habitats that are Polylepis woodlands, unavoidably resulting in some limitations associated with our data. To gather the data to evaluate species richness in Polylepis patches, we used an adapted version of the MacKinnon list method, previously shown to be effective to survey the avifauna of species-poor and remote habitats such as Andean Polylepis forests (poulsen et al. 1997;o'Dea et al. 2004;maCleoD et al. 2006). While this method proved useful to survey birds in our study area, numerous visits of the patches were required to approach the asymptomatic species richness and reduce the uncertainty associated with the use of species richness estimators (maCleoD et al. 2011). Due to the long traveling distances between watersheds, we were also unable to randomize our visits to the different sites, a condition necessary to generate independent visits among sites. Additionally, and due to difficulty of access, we could not survey three of the watersheds in the area, ESP, ICA and LIR, for which we relied on lists generated by a local expert. We found that one of these areas, LIR, appeared to contain distinct communities as compared to the other watersheds (Fig. 3b), although this difference was not explained by any of the environmental factors studied. More data would be needed to determine if this pattern is due to other environmental factors not captured by our own data or to eventual observer biases. The latter seems however unlikely as this pattern was only observed in one of the three areas surveyed by the local expert.
We conclude that Polylepis tree species, habitat quality, topographic complexity and water availability influence avian species communities in Polylepis woodlands. While responses to forest fragmentation likely are species-specific, our results show that higher quality patches characterized by older Polylepis trees, denser canopies and less signs of human activity, are favored by Polylepis-associated bird species. Our results emphasize the importance of local topography in determining the bird species richness that a forest fragment may sustain, through topographical controls on moisture distribution and illumination, and thus on forest productivity and ultimately on the resources available to the avifauna. Conservation and reforestation schemes in similarly fragmented montane forests should therefore account for the topographic characteristics of forest patches, as topography can contribute to habitat quality of patches.
This relationship may become even more important when considering the conservation of such forests and their biodiversity in the face of current climate change trends suggesting that the High Andes may become drier in the future (gareCa et al. 2010;Chevallier et al. 2011;marCora et al. 2013). Furthermore, the establishment of exotic trees close or within native forest fragments should be discouraged. Finally, we suggest that conservation and management efforts should be directed towards both Polylepis forest types (P. lanata and P. subtusalbida fragments) as they contain different bird communities, and that small forest patches also contribute to local bird communities, as habitat quality (i.e. dense forests with tall trees) proved to be more important than fragment size itself, indicating that small fragments should also be conserved.

Supplementary material S1 Detailed technical information: calculation of the TPI
The Topographic Position Index (TPI, Weiss 2001), a measure of local topographic variation, was calculated as the difference between the elevation of the focal pixel and the average elevation of all surrounding pixels within a pre-determined radius (here: 25 m). Pixels with large absolute TPI values (larger than the SD value of the TPI layer, following Weiss 2001) were categorized as being positioned in a depression (negative values) or on a ridge (positive values) while pixels with small absolute TPI values were categorized as belonging to slopes. We used slope values to distinguish between pixels located on a steep slope when their slope angle exceeded 34 degrees, as this corresponds to the maximum inclination where agriculture is practiced within the TNP (sanabria siles et al. 2012) or on a gentle slope if otherwise. The proportion of pixels in a patch belonging to each of these four categories (i.e., depression, ridge, steep slope, gentle slope) was used as indicator for local topographic variation.     S1 Akaike's Information Criterion weights (AIC W ), values of the coefficient and standard errors (SE) of the variables included in the Generalized Linear Mixed Models analysis with Chao2 and ICE total species richness estimates at watershed and patch levels. Variables with AIC W > 0.50 are shown in bold. * indicates variables which have a coefficient value being larger than their standard errors.

TABLE S2
(continued on next five pages). Bird species with their code names observed in Polylepis fragments of the Tunari National Park with their IUCN conservation status (LC = Least Concern, NT = Near-Threatened; EN = Endangered) and the Important and/or Endemic Bird Area (IBA/EBA) to which they belong.       (continued on five next pages). Bird species observed in the Tunari National Park, reported per watershed (in bold) and Polylepis patch ID (second row), when applicable. Complete bird names and codes are reported in Table S2