Tropical mountains are global hotspots for birdlife. However, there is a dearth of baseline avifaunal data along eleva-tional gradients, particularly in Africa, limiting our ability to observe and assess changes over time in tropical montane avian communities. In this study, we undertook a multi-year assessment of understory birds along a 1,750 m elevational gradient (1,430-3,186 m) in an Afrotropical moist evergreen montane forest within Ethiopia's Bale Mountains. Analyzing 6 years of systematic bird-banding data from 5 sites, we describe the patterns of species richness, abundance, community composition, and demographic rates over space and time. We found bimodal patterns in observed and estimated species richness across the elevational gradient (peaking at 1,430 and 2,388 m), although no sites reached asymptotic species richness throughout the study. Species turnover was high across the gradient, though forested sites at mid-elevations resembled each other in species composition. We found significant variation across sites in bird abundance in some of the dietary and habitat guilds. However, we did not find any significant trends in species richness or guild abundances over time. For the majority of analyzed species, capture rates did not change over time and there were no changes in species' mean elevations. Population growth rates, recruitment rates, and apparent survival rates averaged 1.02, 0.52, and 0.51 respectively, and there were no elevational patterns in demographic rates. This study establishes a multi-year baseline for Afrotropical birds along an elevational gradient in an under-studied international biodiversity hotspot. These data will be critical in assessing the long-term responses of tropical montane birdlife to climate change and habitat degradation.

Statistical Analyses

Community-level Analyses

To test whether our survey effort had adequately surveyed the local bird community, we calculated rarified species accumulation curves across sampling days for each site, based on observed and expected (sample-based rarefaction) species richness (Colwell et al. 2012) using the “exact” method of the specaccum function from the R package VEGAN (Oksanen et al. 2019). Since our species accumulation curves did not reach asymptotes for species richness, observed species richness likely does not capture true species richness. We, therefore, used sample-size-based rarefaction and extrapolation (R/E) of Hill numbers (the effective number of species, which integrates species richness and relative abundances; Chao et al. 2014). Sample-size-based rarefaction and extrapolation of Hill numbers is an emerging approach used to standardize and compare estimates of diversity between samples (see Cox et al. 2017, Fair et al. 2018, Baumel et al. 2018, Chao et al. 2019, Debela et al. 2020). Specifically, we used this framework to estimate two values of Hill number 0 (i.e. estimated species richness). First, we calculated standardized species richness. We used the function iNEXT from the R package iNEXT (Hsieh et al. 2016) to calculate R/E curves, standardizing our curve parameters to a maximum of 1,000 individual bird captures (endpoint = 1,000), knots = 500, and a bootstrap replication of 1,000 (nboot = 1,000). From these curves, we provide standardized estimates of species richness based on the sampling of 1,000 individuals at each site. We also estimated asymptotic species richness using the function ChaoRichness from the package iNEXT (Hsieh et al. 2016). Although the asymptotic species richness
is an estimate of true species richness, in practice, reaching an asymptote can take a long time and a lot of sampling. We then plotted the R/E curves of standardized species richness (i.e. over 1,000 individuals) for each site as a function of sample size using the function ggiNEXT (Hsieh et al. 2016). We also visualized asymptotic species richness by setting the endpoint of the iNEXT function to 10,000 individuals. 

Next, we assessed the spatial and temporal patterns in observed species richness and guild-specific captures. For guild-specific captures, we identified the primary diet and habitat association of each species using a global dataset of avian ecological traits (Table 1; see Şekercioğlu et al. 2004, 2019 for a description of the dataset) and summed captures for each separate guild based on either primary diet or habitat. We restricted our analyses to guilds that had ≥40 captures and ≥5 species over the study period and modeled each guild independently. We chose a ≥40 capture threshold as our cutoff between infrequently and frequently encountered species. Most species above this threshold were recorded each year and more than once or twice in each year (the few species that were not recorded each year were recorded multiple times in the other years), whereas individuals under this threshold tended to have few captures across more than one year. We chose a ≥5 species threshold for the guild models to ensure that results for these metrics represented more than a few species.

We constructed models comparing each response variable (observed species richness, dietary, and habitat guild-specific captures) as a function of the site, and
included the number of survey days per site and year (Table 2) as a covariate to control for the variation in the sampling effort. We used generalized linear models (GLMs) for species richness and guild-specific captures, as these represent count data. Within the GLMs, we used a Poisson error structure for species richness, and for guild-specific captures, we used a quasi-Poisson error structure to account for over-dispersion in the count data. To assess changes in the bird community over time, we ran an additional model for each response variable that contained year and site, with a year * site interaction (error structures were applied as above). We tested the significance of the explanatory variables in the GLMs with an analysis of deviance.

We assessed species dissimilarity between sites along the elevational gradient by calculating the Sørenson dissimilarity index (S₈) for pairs of sites adjacent to each other along the elevational gradient, as well as for Chiri-1430 and Dinsho-3186 at either end of the gradient. S₈ can range from complete dissimilarity (S₈ = 1) to complete similarity (S₈ = 0). This dissimilarity can be further decomposed into turnover and nestedness, which we calculated using the function beta.pair in the package betapart (Baselga et al. 2020). Finally, to compare community composition (captures of different species, weighted by abundance), we ran a Principal Coordinate Analysis (PCoA) based on a Bray-Curtis dissimilarity matrix (Legendre and Legendre 2012). A PCoA extracts the greatest orthogonal axes of variation in community composition, plotting them in multidimensional space such that more similar communities are closer to each other in Euclidean space. We extracted the first two axes from the PCoA that represent the greatest variation in community composition.

Species-level Analyses

As a proxy for species abundance (Dulle et al. 2016), we calculated species-specific captures (the number of captured and recaptured individuals of a particular species) per site and year for the most frequently-captured species (≥ 40 captures over the study period). To assess the variation in species’ elevational distributions, we calculated the mean elevation at which each species was detected each year (hereafter “mean elevation”) for frequently-captured species that were detected at least once in every year of the study. Smaller range shifts in tropical birds are more detectable when analyzing mean elevational occurrence rather than the changes in upper or lower range boundaries, as the position of range boundaries is strongly dependent on the sampling effort (Shoo et al. 2006).

We regressed both species-specific captures (in a GLM with a quasi-Poisson error structure) and mean elevation (in a simple linear model) against year. Since the Dinsho-3186 site was located far from the other sites, we decided to re-run the species-level analyses with Dinsho-3186 data removed. The results remained similar with Dinsho-3186 excluded (Supplemental Material Tables S2 and S3) and, therefore, we retained Dinsho-3186 data in the analyses to increase our statistical power. Additionally, we compared our elevational records for banded birds with those reported in the literature for Ethiopia and the Horn of Africa (Ash and Atkins 2009, Dowsett and Dowsett-Lemaire 2015, Rannestad 2016) in order to assess whether any species were detected outside of their recorded elevational distributions. We used an elevational difference of at least 150 m to indicate whether a species had clearly been recorded in our study higher or lower than previously reported in Ethiopia, a distance previously used to signify extralimital records of birds in Ethiopia (Dowsett and Dowsett-Lemaire 2015). A difference of <150 m could result from chance, whereas a difference >150 m is more likely to result from a systematic change in the elevational range.

At the population level, we used Pradel models (Pradel 1996) implemented with the package RMark (Laake and Rexstad 2012) to estimate the rates of apparent survival (φ), recruitment (F), and realized population growth (λ) while controlling for encounter probabilities (p). φ is the rate at which individuals remain in the population; F is the rate at which new individuals join the population via birth or immigration; and λ is the combined effect of survival and recruitment. A population does not change in size when λ = 1, declines when λ <1, and grows when λ >1. These mark-recapture models cannot distinguish movement in and out of a study area (immigration/emigration) from true birth and survival. However, birds living in tropical mountains are known to have small range sizes (Orme et al. 2006), and tropical understory birds of interior forests have limited dispersal ability (Janzen 1967, Moore et al. 2008, Lees and Peres 2009, Visco et al. 2015, Polato et al. 2018, Sheard et al. 2020), suggesting that the differences in recruitment or apparent survival rates of understory birds at our study sites are unlikely to be highly affected by immigration and emigration.

We initially considered all species with >50 captures. For every species we estimated all 4 demographic parameters in constant (time-invariant) Pradel models, extracting the parameter estimates and their 95% confidence intervals (CI). We then censored any species where the models failed to estimate parameters, or where the 95% CI on the parameter estimates were exceedingly large. Then, to assess whether the population growth rates depended on elevation, we modeled λ as a function of the site for each species. Additionally, we ran models where both λ and p were constrained by site, but these models involved many parameters and they were, therefore, suitable only for those species with the largest sample size.

All statistical analyses and graphing were conducted in R (R Core Team 2020; version 4.0.2, 2020-06-22).

For the individual pages of data pertaining to each individual site, the name of the site includes the elevation in m asl for this location. Otherwise, there is no additional information that one needs to know to use this dataset. There are no missing values that would impede someone from re-running our analyses. Please note that there are multiple tabs of data in this file for the dataset. Please see our manuscript for additional information about the study site and the bird surveys.