Community modeling reveals the importance of elevation and land cover in shaping migratory bird abundance in the Andes.

The tropical Andes are characterized by extreme topographic and climatic complexity, which has likely contributed to their outstanding current species diversity, composed of many range-restricted species. However, little is known about how the distribution and abundance of highly mobile organisms, like long-distance migratory birds, varies across different land covers, elevations, and climatic conditions within the Andes. We conducted 1,606 distance-sampling point counts across the Colombian Andes, spanning elevations from 253 to 3,708 m, a range of precipitation regimes and representative land covers. We then employed a novel application of a multispecies hierarchical modeling approach to evaluate how elevation, local land cover, aboveground woody biomass, cloud cover, precipitation, and seasonality in precipitation shape the abundance of the migratory land bird community in the Andes. We detected 1,824 individuals of 29 species of migratory land birds, six of which were considered incidental in our study region. We modeled the abundance of the remaining 23 species, while considering observer and time of day effects on detectability. We found that both elevation and land cover had an overriding influence on the abundance of migratory species across the Andes, with strong evidence for a mid-elevation peak in abundance, and species-specific responses to both variables. As a community, migratory birds had the highest mean abundance in shade coffee plantations, secondary forest, and mature forest. Aboveground woody biomass did not affect the abundance of all species as a group, but a few showed strong responses to this variable. Contrary to predictions of a positive correlation between abundance and precipitation, we found no evidence for community-level responses to precipitation, aside for a weak tendency for birds to select areas with intermediate levels of precipitation. This novel use of a multispecies model sheds new light on the mechanisms shaping the winter distribution of migratory birds and highlights the importance of elevation and land cover types over climatic variables in the context of the Colombian Andes.

Introduction

The tropical Andes are a biodiversity hotspot with exceptional levels of bird diversity (Myers et al. 2000, Kattan and Franco 2004, Fjeldså and Bowie 2012). The complex topography of the Andes mountains, rapid transitions from dry to moist environments and stable year round climatic conditions have led to a remarkable diversity (Fjeldså and Lambin 1999, Cadena and Klicka 2007, Sedano and Burns 2010). Both physical (e.g., deep valleys) and ecological (e.g., rapid shifts from dry to moist conditions) barriers are thought to have historically constrained the dispersion of birds across space, promoting sharply defined distribution limits and speciation (Graham and Silva 2010, Gutiérrez‐Pinto et al. 2012, Valderrama 2014, Winger and Bates 2015, Hazzi et al. 2018). The tropical Andes, as well as other tropical mountains, are further characterized by a rapid species turnover across elevations with species often occupying narrow elevational ranges (Terborgh 1977, Patterson et al. 1998, Jankowski et al. 2009).

In addition to the remarkable diversity of resident species, the tropical Andes are also recognized for their role in providing non‐breeding refuges for a diverse community of migratory birds that breed in North America (Davidson and Fernández 2011). However, our knowledge of how factors such as topography, forest cover at local and landscape scales, and climate limit the distribution of migratory species in the northern Andes remains limited (but see Céspedes and Bayly 2018). Migratory land birds are highly mobile species capable of flying thousands of kilometers without need to rest, enabling them to overcome the most formidable of geographical barriers such as oceans and deserts (Hedenström 2010, DeLuca et al. 2015, Cano et al. 2020). Therefore, unlike many resident birds, high ridges and low valleys associated with the Andes are not expected to physically limit the potential distribution of migrants in northern South America, although they can influence where intraspecific populations with different breeding origins settle (González‐Prieto et al. 2017).

Because of the mobility of migratory land birds, environmental and land cover gradients may play an important role in shaping the distributions and spatial patterns of abundance. Many migratory species are known to track a given environmental niche throughout the year (Gómez et al. 2016, Zurell et al. 2018) and distributions may be shaped by the presence of similar environmental conditions to those experienced on the breeding grounds. Niche tracking has the potential to lead to species‐specific patterns of distribution and abundance but equally, there may be commonalities across species with respect to preferred elevations or climatic conditions. For example, several studies have found a higher abundance of Nearctic‐Neotropical migratory land birds wintering in moister environments (Sherry and Holmes 1996, Smith and Reitsma 2010, Céspedes and Bayly 2018), with moisture in the form of precipitation being directly linked to resource availability in tropical habitats (Studds and Marra 2007). The life cycle of migratory birds is also strongly associated with seasonal environments that provide surplus energy resources at given times of year (McArthur 1959, Somveille and Rodrigues 2018) that migrants might profit from and use to carry out activities such as breeding, molting, or migrating (Leu and Thompson 2002, Bayly and Gómez 2013). Highly stable environments with limited annual variation in available resources and hence predicted to be close to their carrying capacity, might therefore be avoided by migratory birds, as has been shown for breeding distributions (Somveille et al. 2018). Accordingly, we might expect migratory land birds to seek regions of the tropical Andes with stronger seasonality in precipitation during key periods such as late winter when individuals need to start preparing for migratory departure. Given the strong link between precipitation and productivity in the tropics (Wolda 1978, Poulin and Lefebvre 1992, Studds and Marra 2005), this seasonality should lead to variation in carrying capacity that migratory organisms can exploit.

In addition to broadscale environmental conditions, local factors are expected to influence where birds settle across the Andes, including local land cover type and extent of forest cover at broader scales. Earlier studies have indicated that Neotropical migratory land birds use a range of land covers along a forest to open (agricultural landscapes interspersed with trees) gradient during the non‐breeding season, in contrast to the breeding grounds where they primarily use forested landscapes (Petit 1995, Greenberg and Bichier 2000, Endenburg et al. 2019). Some species show higher abundance in mature tropical forests (Céspedes and Bayly 2018, Bulluck et al. 2019), while others occupy live fences in open agricultural areas more than closed native forest (Valdez‐Juárez et al. 2018). The influence of environmental conditions and land cover on overwintering distributions may also be constrained by latitude. In some cases, selective pressures to minimize the distance between the breeding and wintering grounds (Somveille and Rodrigues 2015) may promote greater settlement in the northernmost Andes. Therefore, measuring and accounting for a latitudinal effect on distributions will allow for a more precise identification of how elevation, land cover and climate influence distributions within the latitudinal range of each species.

Previous attempts to study the distribution of migratory birds across the tropical Andes have been focused on single species using either presence‐only data (Colorado et al. 2012) or, in rare cases, abundance data (Céspedes and Bayly 2018). In this study, we identify the factors that drive the distribution and abundance of the community of migrants across the three chains of the Andes in Colombia. The complex topography of the Andes generates strong abiotic and biotic gradients at relatively small spatial scales, making them an ideal system to examine the ecological drivers of migratory bird abundance. We specifically examined the relative influence of elevation, latitude, local land cover type, landscape forest cover, mean precipitation, precipitation seasonality, and mean cloud cover on the abundance of migrants throughout the Colombian Andes. We utilized a novel hierarchical community model (Sollmann et al. 2016) for this analysis. Hierarchical community models have become important tools for the analysis of multispecies data sets with the sharing of information across species but also the estimation of species‐specific relationships (Dorazio and Royle 2005, Gelfand et al. 2005, Zipkin and Dewan 2009). This framework allows for the inclusion of rare species that would otherwise have insufficient information for individual analysis and to identify how the community as a whole responds to variables of interest. The Sollmann et al. (2016) model that we used was developed for distance sampling (Buckland et al. 2001) and originally applied to line transect data; in this study, we adapted it for a point count data set. Through the use of a community model, we were able to evaluate whether migratory species as a group show common responses to geography, land cover, and climatic conditions or if they tend to show species‐specific patterns.

Methods

Study localities

We conducted point count surveys at 1,606 stations located within 62 study sites throughout the Colombian Andes (Fig. 1; for site list see Data S1: Locality_information_AndesMigrants.csv). Study sites were primarily located in the three Colombian Andes ranges (Eastern, Central, and Western Andean cordilleras) with a few additional sites located in the Upper Magdalena region, where the Central and Eastern Andes merge. We selected sites to cover an elevation gradient from 253 to 3,708 m and to ensure inclusion of all representative land covers in the Colombian Andes: forest, live fences, coffee plantations, and paramos (see Point and land cover description for land cover classification), while allowing for constraints associated with ease of access and the availability of trails. We did not sample pastures without trees or bushes because the study species (primarily forest dwelling migratory land birds) do not occur in open pastures (N. J. Bayly, personal observation).

Fig. 1

Distribution of point counts across the Colombian Andes and by elevation and habitat type. (A) Location of counting stations in the Colombian Andes. Black dots represent each of 67 study sites and blue circles vary in size based on the number of point count stations per site (range 1 to 133). Darker gray areas on the map indicate higher elevations. (B) Histogram of number of point count stations (left y‐axis). by elevation. Lines show the density distribution of points for each habitat separately (right y‐axis).

Point count protocol

We conducted 5‐minute passive point counts with distance sampling (Buckland et al. 2001) at each study site between 06:00 and 10:40 from 1 December through 15 March. Distance sampling requires a number of assumptions including (1) all individuals are detected at the survey point, (2) individuals are randomly distributed (though effects of clustering can be modeled), and (3) distances to individuals are accurately measured. We defined survey dates to maximize the inclusion of overwintering birds and exclude passage migrants during the migration periods. To establish point count locations, we paced out 100 m or more between points, which, given that most points were located on steep Andean slopes, resulted in horizontal distances of around 80 m when measured by GPS but distances over the ground of 100 m or greater. We conducted surveys during the winters of 2013, 2015, 2016, 2018, and 2019, however, 96% of surveys were carried out over three winters (2015, 2016, 2019). The number of point count stations per site was highly variable (range = 1–133, mean = 25.9; Data S1) and largely depended on the size of a site (i.e., patch size of forest or coffee plantations). Six sites were sampled in multiple wintering periods and in each year point locations were different, but not necessarily separated by 100 m over the ground among years. We did not conduct counts if moderate to heavy rain or dense mist was present.

During each point count, we recorded all Nearctic‐Neotropical migratory land birds seen or heard within a radius of 40 m, although the majority of individuals were detected within the first 25 m (86%). For each individual detected, we estimated the horizontal distance from the observer (i.e., from the center of the point) to the nearest meter by eye, occasionally verifying by pacing or using a range finder. Point counts were carried out by four observers: Nicholas Bayly, Laura Céspedes, Camila Gómez, and Natalia Cano. All observers were trained in the visual and aural (using calls) identification of migrants that overwinter in the Colombian Andes, as well as in the estimation of horizontal distance, by comparing among observers and through calibration using a range finder.

Point and land cover description

At each point count station, we recorded coordinates and elevation using a GPS with a precision of 10 m or higher. We categorized land cover at each point according to the following classes (see Data S1): mature forest (938 points), secondary forest (322), forest fragments (64), shade coffee (192), sun coffee (29), live fences (32), and paramos (29). Mature forest was defined as having a canopy height >15 m but above 2,250 m the assignment was more subjective, as elfin forests were common and mature stands rarely exceeded 15 m in height. For example, we classified points as mature forest at high elevations when trees with estimated diameter at breast height (DBH) >30 cm were present and there was a high abundance of bromeliads (Dunn 2000). We used the category forest fragment for points located in small patches (<1 ha) of forest with high canopy height (>15 m); these patches were usually associated with small streams or rivers. To classify coffee plantations, we followed Moguel and Toledo (1999) and Perfecto et al. (2005) and defined polycultures and rustic plantations with a canopy cover >30% as shade coffee and monocultures with canopy cover <30% as sun coffee. We estimated canopy cover by looking through inverted binoculars and estimating the visual field that was covered by vegetation. Live fences corresponded to points in linear vegetation in open areas, usually dominated by cattle pastures with occasional shade trees. We classified points above 3,300 m with grassy vegetation, dispersed bushes, and Espeletia plants as paramo (high elevation grasslands characteristic of the northern Andes). It is noteworthy that, except for the paramo category, which naturally lacks trees, even open land covers such as live fences and sun coffee usually occurred in a matrix that included small patches of forest or secondary growth.

Acquisition of environmental data

To extract values of precipitation (monthly; Fick and Hijmans 2017), cloud cover (monthly; Wilson and Jetz 2016), and aboveground live woody biomass density (Baccini 2012) for each point count station, we used/created raster layers with a 1‐km resolution (woody biomass density averaged from 30‐m resolution to get a 1‐km resolution; Appendix S1: Fig. S1). The aboveground live woody biomass layer describes the amount of biomass corresponding to trees per unit of area and here we used it as a proxy for forest cover and maturity (i.e., areas with continuous mature forest should have the highest values). We constructed the layer using images for the year 2000, meaning that biomass values when point counts were conducted might have been slightly different due to deforestation or restoration in recent years. However, variation in forest cover indicated that the rate of forest loss was low in the study area between 2000 and 2019 (around 3.4% for elevations between 1,000 and 2,300; N. J. Bayly, unpublished data).

In contrast to previous studies, we used precipitation and cloud cover to represent productivity rather than NDVI, in part due to the inherent error in estimating NDVI in cloudy environments (e.g., Andean cloud forests) and also because of a growing body of work that links winter precipitation/moisture to body condition, overwinter survival, migratory condition, and subsequent reproductive success in migratory land birds (Norris et al. 2004, Studds and Marra 2007, González et al. 2020). To describe average values for the non‐breeding period (November–April), we generated new raster layers for precipitation and cloud cover (note annual and non‐breeding values were correlated; data not shown). For precipitation we also calculated a seasonality variable to describe the degree to which rainfall increased toward the end of the wintering period, immediately prior to spring migration when the energy requirements of migratory birds are elevated: Seasonality = ((mm of precipitation in March and April) – (mm of precipitation in December through February))/(mm of precipitation in December through April). As calculated here, values for seasonality will be higher for areas with greater amounts of precipitation at the end of the wintering season (March and April), and lower for areas with greater amounts of precipitation at the beginning of the wintering season (December through February). The denominator was included to account for absolute differences across sites in total precipitation through the winter period. We manipulated raster layers using functions implemented in the raster (Hijmans and Etten 2012), rgdal (Bivand and Keitt 2017), and gdalUtils (Greenberg and Mattiuzzi 2018) packages for R (R Development Core Team 2019). Correlation between environmental variables included in the model was relatively weak overall (Appendix S1: Fig. S2).

Hierarchical community distance‐sampling model

We used hierarchical distance sampling (Buckland et al. 2001) in a community modeling approach to evaluate the variables influencing the distribution and abundance of the migrant community (code in Data S2, data used as input in Data S3). Here we use the term community to refer to the set of species that winter across the Andes following the definition provided in Ricklefs (2008). We excluded some species from analyses that were detected in very low numbers in our point counts (Appendix S1: Table S1) and whose wintering distribution is primarily in the lowlands of the Caribbean, Orinoquia, or the Amazon rather than Andean (e.g., Baltimore Oriole Icterus galbula, Blackpoll Warbler Setophaga striata).

With distance sampling, we assume perfect detection for all individuals at the survey point (i.e., radius = 0 m) and a declining detectability with distance defined by a function g. We used a half‐normal detection function with x being the distance from the center of the radius and σ being the scale parameter describing the shape of the function. Individual detections were grouped into 5‐m distance bands from 0 m up to and including 40 m.

The detection parameter for each species s and station j can be modeled as a function of covariates (e.g., Marques and Buckland 2003) where α0 is the detection intercept for species s and is a vector of species‐specific coefficients describing the relationship between detection probability and the observation covariates Y at station j. These covariates included observer effects, time of day, and canopy height. Note that we were unable to include species‐specific observer effects.

The observed number of individuals n of species s at each station j is linked to the true abundance N via the detection probability p.t

The abundance parameter, N, was assumed to follow a Poisson distribution (Royle and Dawson 2004) with expected value λ defined by the covariates describing elevation, land cover, latitude, and environmental conditions where β0 is the mean abundance (i.e., intercept) for species s and is a vector of coefficients for the relationship of species s to the set of covariates X measured at each point count station j.

For a hierarchical model with s = 1, 2, …, M species, the species‐specific parameters are described by hyperparameters representing the community where μβ0 and σ2 β0, respectively, describe the mean and variance of the community response to each of the k = 1 – n variables described by the β parameters. This approach is similar to including a normally distributed random effect for each species in the model (Sollmann et al. 2016).

We initially considered four variables that might influence detection probability: observer (categorical with four levels), linear time of day (linear term only), curvilinear time of day (linear and quadratic terms), and canopy height. We included the time‐of‐day variables to test whether (1) detection increased or decreased consistently through the morning (linear) or (2) whether detection was highest or lowest mid‐morning (curvilinear). We included canopy height as an index of forest cover to incorporate variation in detectability among land cover types. Time of day and canopy height were allowed to vary by species, but we could only include observer effects as a common response across species due to poor model convergence in the observer by species interaction model. We suspect that this poor performance was because the less common species were often observed by only a single observer. After a preliminary analysis with all detection and abundance covariates, we found no evidence for a curvilinear time‐of‐day effect for any species and therefore, the final set of models included only effects of observer, time of day (linear), and canopy height on detectability.

The full set of abundance covariates included elevation, latitude, woody biomass at a 1‐km scale, land cover class at the point count station, mean precipitation, seasonality in precipitation, and mean cloud cover. For woody biomass and environmental variables, we considered a linear model (linear term) and a curvilinear model (linear and quadratic terms) because it is possible that species select intermediate values. For elevation, we only included the curvilinear model to best identify the elevation at which community and species‐specific abundance peaks. We standardized all continuous covariates to a mean of 0 and a standard deviation of 1 prior to analysis. For the land cover classes, we set mature forest as the basal group and thus, all coefficients reflect the difference in the expected log abundance of each species and the community mean in each land cover class relative to that in mature forest.

We implemented the hierarchical community model with a Bayesian framework using JAGS (Plummer 2003) in R Version 3.6.0 (R Development Core Team 2019). We used vague priors for all detection and abundance parameters (see code in Data S1). We ran three Markov chains in parallel with inference based on 5,000 iterations after an initial 1,000 adaptive iterations and a burn‐in of 14,000 iterations. The final sample of 5,000 iterations was thinned by five and we report the posterior means and the lower and upper bounds of the 90% interval around parameter estimates. We assessed model convergence using the R‐hat statistic, which is a measure of the variance within vs. among Markov chains; convergence was assumed when R‐hat values were <1.1 (Gelman et al. 2004).

We examined model fit using Bayesian P values for the abundance and detection components following Sollmann et al. (2016). For the abundance P value, we generated new data N for every species and station given the expected values of λ. The Bayesian P value was then based on a comparison of the observed residuals calculated between the data N and λ vs. the simulated residuals between the data N and λ. The fit of the detection component is similar with a comparison of the observed number of individuals y for the distance classes. In this case, we simulated a new set of distance classes and compared the residuals for the observed vs. expected detections given the half‐normal distance function. For well fit data, the observed residuals should be no more extreme than the simulated residuals resulting in a Bayesian P value near 0.5. Values that are less than 0.1 or greater than 0.9 suggest lack of fit.

Estimation of species and community density

We used the model covariates to estimate densities for all species and the full community in each land cover class. For all classes except paramo, we based estimates on the average predicted abundance between 1,000 m and 2,000 m elevation using the coefficients for elevational effects across that gradient. This elevational range was chosen to predict densities because all land cover types except paramo commonly occur within this range. For paramo, we predicted densities using elevations >3,275 m, the lowest elevation at which paramo occurred across study sites. We scaled densities for mature forest, second growth forest, shade coffee, sun coffee, and paramo to the expected number of individuals per square kilometer of that class. We used a different approach to estimate densities for forest fragments and live fences, which do not form continuous habitat blocks but are instead fragments of potential habitat within otherwise open and mostly unsuitable habitat for our study species. Therefore, for these two land cover classes, we estimated their average extent within 1 km2 of open land, based on typical configurations within our study region, and scaled density estimates to the expected extent. While this approach avoids overestimating densities in landscape mosaics, it is prone to inflating densities within strips of vegetation through the application of the distance correction (Hutto 2016). As described previously, forest fragments were typically riparian strips and in Colombia, landowners are required by law to maintain 30 m of forest cover either side of water courses. We therefore set riparian strip width at 60 m. Then, using satellite imagery, we estimated that, on average, there are approximately two riparian strips, each 1 km in length, per 1 km2 within much of the Andean region. This resulted in an estimated extent of 0.12 km2 riparian strip/km2. We used a similar approach with live fences, estimating an average width of 10 m and ~5 1‐km fences/km2, which gives an average extent of 0.05 km2 live fence/km2. We acknowledge that some suitable habitat may ordinarily be available in these landscapes surrounding fragments and live fences and therefore the estimates in these two land cover classes should be considered minimum estimates.

Results

We detected 1,824 individuals of 29 Nearctic‐Neotropical migrant species across 1,606 point counts. As described earlier, we removed six species from the analysis because they are typically lowland species and were considered incidental observations (Appendix S1: Table S1). This resulted in a final data set of 1,815 individuals of 23 species. The number of detections by species ranged from one (e.g., Alder Flycatcher Empidonax alnorum, Golden‐winged Warbler Vermivora chrysoptera, Scarlet Tanager Piranga olivacea) to 502 for Blackburnian Warbler (Setophaga fusca) (Table 1, Appendix S1: Table S1). After Blackburnian Warbler, the next five most abundant species were Swainson’s Thrush Catharus ustulatus (418 detections), Canada Warbler Cardellina canadensis (273), Acadian Flycatcher Empidonax virescens (204), Summer Tanager Piranga rubra (89), and American Redstart Setophaga ruticilla (80). While we detected 29 species in total across all surveys, including the 6 considered to be incidental; 16 species represented more than 99% of the observations.

Table 2

Posterior summaries for community parameter median coefficients (on log scale) from a community distance‐sampling model for Neotropical migrant birds in the Colombian Andes.

Parameter	Median	Lower 5%	Upper 95%
Elevation	−1.36	−2.02	−0.77
Elevation2	−0.91	−1.23	−0.61
Latitude	0.13	−0.14	0.38
Second growth	−0.08	−0.46	0.27
Shade coffee	0.29	−0.25	0.83
Sun coffee	−0.42	−1.35	0.43
Forest fragment	−0.08	−0.93	0.78
Live fence	−1.04	−2.84	0.34
Paramo	−7.35	−20.52	1.32
Tree cover	−0.05	−0.28	0.22
Mean precipitation	−0.09	−0.33	0.14
Precipitation seasonality	−0.02	−0.23	0.17
Cloud cover	−0.02	−0.17	0.14

Lower 5% and Upper 95% are the lower and upper bounds on the 90% credible interval. Mature forest was the basal categorical group and therefore categorical land cover responses are effects relative to mature forest. See Appendix S1: Table S3 for coefficients from a model with quadratic terms for tree cover and climate variables.

Variation in abundance across the elevation gradient

Elevation was the strongest predictor of relative abundance (hereafter “abundance”), explaining a significant amount of variation at both the community level and the individual species level (Table 2, Fig. 2). As a group, migratory birds showed a common response to elevation, leading to a mid‐elevation peak in predicted abundance of 1,847 individuals/km2 (90% CI 1,663, 2,045) based on estimates in mature forest at elevations of approximately 1,470 m (standardized log scale = −1.36 [90% CI −2.02, −0.77], = −0.91 [−1.23, −0.61], Table 2, Fig. 2A). Abundance declined on either side of this peak reaching a low of 328 individuals/km2 (90% CI 162, 617) at 295 m, just above the minimum surveyed elevation (253 m), and 21 individuals/km2 (90% CI 3, 62) at 3,580 m, near the highest elevation surveyed (3,610 m). While community abundance showed a clear mid‐elevation peak (Fig. 2A), individual species responses to elevation did vary (Fig. 2B, C). For example, Bay‐breasted Warbler Setophaga castanea peaked at around 480 m, well below the community mean, while Blackburnian Warbler and Black‐and‐white Warbler Mniotilta varia both peaked above the community mean. It is noteworthy that the true variation in abundance in relation to elevation might deviate from a quadratic curve. For example, curves might have a mid‐elevation peak but might not be symmetric on either side of the peaks, as tends to be predicted from a quadratic model. Nevertheless, we found that raw data (i.e., proportion of detections per elevation band) followed a similar trend to model predictions (Appendix S1: Fig. S3).

Table 1

Posterior summaries of the estimated densities (individuals/km2) of Neotropical migrants in the Colombian Andes.

Species and no. detections	Mature forest (1024)	Second growth (331)	Shade coffee (199)	Forest fragment (64)	Sun coffee (31)	Live fence (32)
Blackburnian Warbler (502)	369 (305, 439)	454 (343, 568)	797 (615, 998)	34 (18, 55)	469 (175, 859)	33 (16, 56)
Swainson’s Thrush (418)	378 (312, 447)	495 (377, 638)	280 (211, 356)	23 (10, 37)	112 (22, 239)	1 (0, 4)
Canada Warbler (273)	400 (317, 483)	218 (135, 315)	207 (135, 285)	10 (2, 23)	134 (27, 290)	1 (0, 4)
Acadian Flycatcher (204)	148 (112, 188)	183 (117, 257)	97 (61, 140)	15 (6, 27)	47 (2, 121)	2 (0, 7)
Summer Tanager (89)	52 (35, 74)	45 (22, 74)	44 (21, 70)	4 (1, 8)	48 (8, 120)	2 (0, 5)
American Redstart (80)	108 (68, 155)	93 (43, 151)	64 (29, 109)	6 (1, 16)	32 (0, 111)	2 (0, 6)
Tennessee Warbler (64)	16 (6, 29)	23 (6, 49)	115 (58, 195)	2 (0, 7)	16 (1, 62)	1 (0, 3)
Mourning Warbler (42)	16 (7, 30)	21 (7, 45)	60 (24, 106)	7 (1, 17)	29 (2, 109)	6 (1, 15)
Bay‐breasted Warbler (31)	18 (6, 34)	13 (3, 29)	25 (6, 55)	1 (0, 3)	7 (0, 33)	0 (0, 1)
Black‐and‐white Warbler (28)	20 (9, 36)	20 (6, 43)	25 (6, 56)	2 (0, 8)	11 (0, 52)	3 (0, 13)
Cerulean Warbler (24)	24 (11, 45)	17 (4, 38)	46 (15, 92)	1 (0, 6)	11 (0, 52)	0 (0, 2)
Yellow Warbler (16)	4 (0, 11)	3 (0, 10)	9 (1, 26)	2 (0, 6)	5 (0, 26)	1 (0, 4)
Gray‐cheeked Thrush (12)	4 (1, 10)	4 (1, 12)	3 (0, 9)	5 (0, 14)	3 (0, 18)	0 (0, 2)
Western Wood Pewee (11)	4 (1, 9)	4 (1, 11)	11 (2, 25)	0 (0, 2)	11 (0, 60)	0 (0, 1)
Olive‐sided Flycatcher (6)	3 (0, 6)	3 (0, 8)	2 (0, 8)	0 (0, 1)	4 (0, 25)	0 (0, 0)
Northern Waterthrush (5)	2 (0, 6)	2 (0, 6)	3 (0, 12)	1 (0, 5)	1 (0, 11)	0 (0, 1)
Yellow‐throated Vireo (3)	1 (0, 3)	1 (0, 3)	2 (0, 9)	0 (0, 1)	1 (0, 5)	0 (0, 0)
Rose‐breasted Grosbeak (2)	1 (0, 2)	0 (0, 2)	2 (0, 6)	0 (0, 1)	0 (0, 3)	0 (0, 0)
Eastern Wood Pewee (1)	0 (0, 2)	0 (0, 2)	1 (0, 3)	0 (0, 1)	0 (0, 3)	0 (0, 0)
Alder Flycatcher (1)	0 (0, 2)	0 (0, 2)	0 (0, 3)	0 (0, 0)	0 (0, 3)	0 (0, 0)
Golden‐winged Warbler (1)	0 (0, 2)	0 (0, 2)	0 (0, 3)	0 (0, 2)	0 (0, 3)	0 (0, 0)
Scarlet Tanager (1)	0 (0, 2)	0 (0, 2)	1 (0, 4)	0 (0, 1)	0 (0, 3)	0 (0, 0)

Densities were estimated for each habitat and are the average for each species across elevations ranging from 1,000 to 2,000 m with their respective 90% confidence intervals in brackets. This range was chosen to predict densities across the elevational range at which all habitat types commonly occur. No individuals of these species were detected in paramo. Species are ordered by the number of detections (see value in parentheses after species name). Values under habitat names refer to the number of point counts by habitat type.

Fig. 2

Estimated abundance in relation to elevation in the Colombian Andes for (A) the community of 23 boreal migrant passerines and (B, C) the 14 most common boreal migrants representing 98% of community abundance. Solid line in panel A shows the median community abundance with dashed lines representing the 90% credible intervals. Individual species abundances were standardized for panels B and C to show them on the same scale. Species four‐letter codes in the legend are as follows: BLBW, Blackburnian Warbler; SWTH, Swainson’s Thrush; CAWA, Canada Warbler; ACFL, Acadian Flycatcher; AMRE, American Redstart; TEWA, Tennessee Warbler; SUTA, Summer Tanager; MOWA, Mourning Warbler; BAWW, Black‐and‐white Warbler; BBWA, Bay‐breasted Warbler; CEWA, Cerulean Warbler; YWAR, Yellow Warbler; GCTH, Gray‐cheeked Thrush; WWPE, Western Wood‐Pewee.

Influence of local land cover and landscape woody biomass on abundance

Both community abundance and species‐level abundance varied strongly in response to land cover class. Based on estimates between 1,000 m and 2,000 m elevation, community abundance was highest in shade coffee (1831 individuals/km2, 90% CI 1,576, 2,100), followed by secondary forest (1,633, 90% CI 1,402, 1,774), mature forest (1,587, 90% CI 1,442, 1,736), sun coffee (1,046, 90% CI 614, 1,531), forest fragments (124, 90% CI 92, 162), and live fences (58, 90% CI 34, 85). Estimates of abundance for paramos above 3,300 m were the lowest at 0.07 individuals/km2 (90% CI 0, 40). There was considerable variation across species in abundance by land cover type (Table 1, Appendix S1: Table S2). While some, such as Blackburnian Warbler and Summer Tanager, were common in both forested and semi‐open land covers, others were either more abundant in closed mature and secondary forest habitat (e.g., Canada Warbler, Swainson’s Thrush), or regularly used semi‐open land covers (e.g., Mourning Warbler Geothlypis philadelphia, Yellow Warbler Setophaga petechia; Table 1, Appendix S1: Table S2).

Community mean abundance was not influenced by the amount of aboveground live woody biomass (i.e., forest cover, hereafter “woody biomass”) at the 1‐km2 scale (Table 2) but several species displayed individual responses (Fig. 3). American Redstart, Bay‐breasted Warbler, Summer Tanager, and Yellow Warbler all showed negative responses to landscape‐scale woody biomass while Gray‐cheeked Thrush Catharus minimus showed a near‐significant positive relationship, with the lower 5% bound only slightly overlapping 0. Swainson’s Thrush showed a curvilinear relationship with abundance peaking in landscapes with intermediate woody biomass (Appendix S1: Fig. S4). The average proportion of forest cover across all our study sites was approximately 75%, with more than 50% forest cover in >95% of sites, at least for the year 2000 (Hansen 2013). Thus, we interpret our results as showing that species whose abundance declined with increasing forest cover at the 1‐km2 scale tended to favor semi‐forested sites.

Fig. 3

Coefficients and 95% credible intervals for species relationships with tree cover, mean precipitation, seasonality in precipitation, and mean cloud cover. Species are arranged from positive to negative responses. Species in red and blue are those that showed a significant positive or negative linear relationship, respectively, with the variable based on the lower 5% bound or upper 95% bound excluding 0 for the estimated coefficient (for species in black credible intervals overlapped 0). All results shown are based on a model with a linear relationship but those species with an asterisk showed evidence for a curvilinear relationship based on the coefficient for the quadratic term excluding 0 (see Appendix S1: Table S3 for quadratic model coefficients).

Responses to latitude and climate variables

Latitude was a strong predictor of abundance for some species with abundance increasing to the north for four species (Acadian Flycatcher, Bay‐breasted Warbler, Mourning Warbler, and Tennessee Warbler) and to the south for one species, Canada Warbler (Fig. 3). Contrary to our predictions, we found no community‐level linear relationships with mean winter precipitation, seasonality in precipitation or mean cloud cover across sites (Table 2), however, the community tended to select areas of intermediate winter precipitation indicated by a quadratic term with more than 85% of posterior mass below 0 ( = −0.05 (90% CI −0.34, 0.19), = −0.10 (−0.25, 0.03), Appendix S1: Table S3). Swainson’s Thrush and Acadian Flycatcher showed a preference for intermediate precipitation levels, while Canada Warbler showed a preference for intermediate cloud cover (all 90% CIs for quadratic winter precipitation term were negative, Appendix S1: Fig. S4). Blackburnian Warblers were more abundant in sites with greater mean precipitation, while American Redstarts showed the opposite pattern (Fig. 3).

Model fit and detectability

Bayesian P values indicated that the model fitted the data well for both the detection (P = 0.68) and abundance (P = 0.38) components. The median community detection estimate on the log scale was 2.17 (90% CI 1.96, 2.37) and the back‐transformed estimate indicated that on average species were detected at 8.76 m (90% CI 7.10, 10.70). The mean horizontal detection distance across species ranged from a low of 6.93 m (90% CI 5.09, 9.40) for Yellow Warbler to a high of 12.12 m (9.84, 14.87) for Summer Tanager (based on predictions for Observer 1, Appendix S1: Table S5). We found variation in detectability among observers and with time of day and canopy height for some species (Appendix S1: Table S4). Detectability decreased over the morning for two species, Acadian Flycatcher and Swainson’s Thrush, and increased with canopy height for five species, Acadian Flycatcher, Blackburnian Warbler, Canada Warbler, Gray‐cheeked Thrush, and Swainson’s Thrush (Appendix S1: Table S4). All species that showed a positive relationship between detectability and canopy height selected forested land covers and we suspect that this response was due to the greater openness of the understory and subcanopy in taller forests.

Discussion

The Colombian Andes are both a global biodiversity hot spot (Pimm et al. 2014) and an important overwintering region for a community of Nearctic‐Neotropical migrant birds (Davidson et al. 2011), many of which are in decline (Rosenberg 2019, Wilson et al. 2019). Using a recently developed hierarchical distance‐sampling approach (Sollmann et al. 2016), we provide one of the first multispecies‐level assessments on how geographic, environmental, and land cover variables combine to determine the distribution of the migrant community across the complex topography of the Colombian Andes. Our results highlight the strong influence of elevation in shaping migrant distributions, with the community as a whole showing a mid‐elevation peak in species abundance between 1,300 m and 1,700 m. While the community concentrated at mid elevations, we observed considerable local habitat segregation among species across natural and agricultural land covers varying in degree of openness. Species also varied in preferences to tree cover and local climatic conditions with the result that the community as a whole showed no particular relationship to these variables, contrary to our expectation of a positive linear relationship with moisture, represented by our precipitation and cloud cover variables, or with seasonality. However, when we included quadratic terms in our model, there was a tendency for migratory species to be more abundant in areas with intermediate levels of winter precipitation within the Andes.

Variation in abundance with elevation

Resident bird species that occur in the Andes are often specialized and restricted to narrow elevational bands resulting in high community turnover across the elevation gradient (Terborgh 1977, Patterson et al. 1998, Jankowski et al. 2009). Here, we show that long‐distance migrant land birds fit this general pattern and tend to occur within species‐specific and restricted elevation ranges across the Andes during winter (Appendix S1: Fig. S3). In contrast to many resident species, Neotropical migrants are highly mobile and can potentially tolerate a variety of different environmental conditions and land cover types across the annual cycle (Petit et al. 1995, Zuckerberg et al. 2016) and it is unclear why we observed clear elevational patterns in migratory land birds. These patterns might arise because migrants track a specific climatic niche across the year, such that climatic similarities between temperate breeding grounds and specific elevation ranges in the Andes result in the observed patterns. Evidence for year‐round niche tracking exists for long‐distance migrant birds, including several that represent the community in this study, particularly within the Parulidae family (Nakagawa et al. 2004, Gómez et al. 2016). However, other long‐distance migrants are more consistent with niche switchers (Nakagawa et al. 2004, Dufour et al. 2020). The apparent elevational zonation of migratory birds in the Andes might also be explained by biotic interactions, as has been proposed for resident species in mountain systems (Jankowski et al. 2013), especially competition (Jankowski and Robinson 2010, Freeman and Class Freeman 2016).

Migratory birds, as a group, showed an apparent mid‐elevation peak in richness as is the case for resident species on some Andean slopes and in many other mountain systems (Rahbek 1995, 1997, McCain 2009, Price 2014; but see Quintero and Jetz 2018, Pigot et al. 2016). Mid‐elevations may represent intermediate zones where the distribution of multiple species tend to overlap more (Colwell and Lees 2000, McCain and Grytnes 2010). For example, species like the Blackburnian Warbler and Bay‐breasted Warbler, peaked in abundance at higher and lower elevations relative to the community average, yet the distributions of both overlapped at mid elevations occupied by most of the species. We acknowledge, however, that the number of species studied here is small and that analyses were not specifically focused on richness, but on abundance. The higher total abundance of the migrant community at mid elevations may reflect a higher carrying capacity whereby mid‐elevation habitats can support more individuals per unit area during the non‐breeding season (Quaintenne et al. 2011). A plausible explanation for a higher mid‐elevation carrying capacity is a greater abundance of insects, which appears to be driving the mid‐elevation peak in abundance and diversity of birds, and particularly insectivores, in the eastern Himalayas, for example (Price et al. 2014, Schumm et al. 2020). To better understand what drives this mid‐elevation peak in abundance, future studies should also consider abundance patterns in resident species that potentially compete for resources with the migratory species discussed herein.

Abundance by land cover types

Although the distribution and abundance of multiple Neotropical migratory land bird species overlapped at mid elevations, there were marked differences across the community in the use of local cover types. Some species were mostly restricted to forested land covers (e.g., Canada Warbler, Swainson’s Thrush), others primarily used semi‐open agricultural land covers (Mourning Warbler, Yellow Warbler) and some, such as Blackburnian Warbler, were more generalist and common across most land cover types (Appendix S1: Table S2). We also found that only one species responded positively to woody biomass in the broader landscape, while several preferred areas with intermediate values and a few others declined in abundance as woody biomass increased. It is important to note, however, that we focused our sampling effort in relatively forested areas and thus the latter case still represents species preferring semi‐open rather than open landscapes. Nevertheless, our results support earlier research elsewhere in the Neotropics indicating that as a group, Neotropical migrants are often more common in semi‐open, secondary, and edge land cover types, as opposed to closed mature forest (Hutto 1989, Petit et al. 1995, Bakermans et al. 2009). This apparent tolerance to lower forest cover for Neotropical migrants on the nonbreeding grounds differs from the pattern on the breeding grounds, where they favor heavily forested landscapes to a greater extent than residents or short‐distance migrants and are often the first to disappear from the avian community as forest cover is lost due to anthropogenic activities (Endenburg et al. 2019).

It is noteworthy that relative community abundance was similar among shade coffee and native forest (secondary and mature). The high abundance in shade coffee may reflect the fact that quality shade coffee offers a mixture of conditions selected by species that favor forest and semi‐open land covers. Other studies in the Neotropics have also found a similar or higher abundance of migratory land birds in shade coffee than in native forest (Tejeda‐Cruz and Sutherland 2004, Bakermans et al. 2009). We note, however, that shade coffee is grown within a narrow elevational range that overlaps with that preferred by most Neotropical migrants. In contrast, our mature and secondary forest points covered a broader range of elevations above and below this mid‐elevation peak; this greater range may have resulted in a lower average density estimate for mature and secondary forest. It would be valuable to compare community and species densities among these land covers using a reduced set of points within the same elevational and geographic range. This approach could also include interactions between elevation and land cover. For example, the structure of shade coffee may be less variable across an elevational gradient compared to native forest and therefore it would be interesting to test if there is a more limited change in abundance across elevation within shade coffee or other agroforestry types compared to native land covers. While shade coffee supported a high community abundance, some individual species abundances were higher in native mature and second growth forest, including threatened species like the Canada Warbler, whose long‐term population declines are believed to be due to the loss of forest cover on the Andean wintering grounds (Wilson et al. 2018). Shade coffee can provide an important habitat to individuals in the absence of native forest (Céspedes and Bayly 2018, González et al. 2021) but the protection and restoration of native forest needs to be a part of conservation strategies for migratory species as a whole.

Climatic variables and abundance patterns

Macroecological explanations for the seasonal distribution of migratory birds (McArthur 1959, Somveille et al. 2015) have invoked the energy‐efficiency hypothesis, which has strong support from correlative studies and mechanistic models. This hypothesis posits that organisms should be distributed so as to minimize their overall energy expenditure while targeting areas with maximum energy available given the distribution of other species (Somveille et al. 2018). Indeed, the richness of migratory birds on their breeding grounds is linked to the degree of seasonality in resources, as measured through NDVI, and to the carrying capacity of the same environments during the non‐breeding season, reflected in winter harshness (Somveille et al. 2018). Combined, these variables describe regions with the highest variation in carrying capacity and are a proxy for the surplus resources that migratory organisms can take advantage of during the breeding season (Newton and Dale 1996, Dalby et al. 2014). During the non‐breeding season, the large‐scale processes driving distributions are not as clear but are related to minimum temperatures, productivity and distance to breeding grounds (Somveille et al. 2018). In the context of the Andes, productivity might therefore be expected to drive the main settlement patterns in migratory birds.

We expected abundance to be correlated in some way with precipitation or cloud cover, our proxies for productivity, given the direct link between precipitation and insect and fruit abundance in the tropics (Boyle and Shogren 2020). In partial support of this, we found that migratory birds, as a group, tend to be more abundant in areas with intermediate levels of wintering precipitation. It is notable that a significant positive correlation between abundance and precipitation was found in just one species (Blackburnian Warbler), suggesting that the pattern of increasing habitat quality with increasing humidity found in other studies in the Neotropics (Norris et al. 2004, Studds and Marra 2007, González et al. 2020) is not true of the wider community. Indeed, for the American Redstart, which survives and maintains higher body condition in moister habitats on Jamaica (Marra et al. 2015), we found a significant negative relationship between abundance and winter precipitation. Note, however, that we did not sample dry forests and therefore the land cover types included in our study are not comparable to the driest land covers occupied by redstarts in Jamaica. Avoidance of extremely wet regions of the Andes may explain these results and supports recent hypotheses regarding the negative effects of extreme rainfall levels on tropical endotherms (Boyle et al. 2020).

Model caveats

We acknowledge four important caveats in our results. First, the use of a multispecies model allowed us to include species that would otherwise not have had a sufficient number of observations for analysis and this provided strength to our understanding of how the group as a whole responded to gradients of elevation, land cover type, aboveground woody biomass, and climate. However, the hierarchical nature of the model also means that estimates for species with limited data are influenced by the mean community response (Link 1999, Zipkin et al. 2009). A few species in our analysis were very rare with only one or two observations. In these cases, we do not provide any interpretation of their estimated relationships with the variables because those would primarily reflect the community mean. However, the model does appear well able to estimate individual species responses with even a modest number of observations. For instance, we were able to identify the distinctly lower elevational preference for Bay‐breasted Warbler (n = 31 observations). We were also able to identify the positive association of Gray‐cheeked Thrush (n = 12 observations) with aboveground woody biomass (i.e., landscape forest cover), the only species to show this association and reflecting what has been found previously for this species at stopover and wintering sites in Colombia (Bayly et al. 2016, González‐Prieto 2018). Yet, hierarchical community models are novel and additional research on relationships between community and individual species estimates would be valuable areas of further work.

A second related caveat is that the lack of community‐level responses to some variables were driven by the variation in response across species. While this observation reflects the characteristic of the full community of interest in this study, finer level groupings would highlight significant community responses within those groups (Pacifici et al. 2014). For example, based on this study and previous work (Céspedes and Bayly 2018), it would be possible to define groups of species differing in preference to landscape‐scale forest cover and these groups could be used as indicators for efforts such as forest restoration and protection (e.g., Schuster and Arcese 2013).

A third caveat in our work is that even minor biases in the estimate of distance to individuals observed will influence estimates of abundance when extrapolated from the 40 m radius point to the square‐kilometer scale, as will the movement of birds into the counting area during the count period (Raman 2003, Buckland and Marsden 2008). Nonetheless, observers in our study placed considerable emphasis on calibrating distance and ensuring that it was similar among observers, and comparable density and home range estimates suggest our estimates are reasonable. For example, Céspedes and Bayly (2018) estimated densities for Canada Warblers in the Andes up to 400 individuals/km2, while home range size for Swainson’s Thrush in Andean forest was estimated at 0.3 ha (McElaney 2019), which is equivalent to a density of 333 individuals/km2 assuming no overlap in home range. However, we still recommend that the density estimates in this study be considered relative and not absolute estimates (e.g., among land cover types, Table 1).

Finally, we acknowledge that species abundances may vary among years and that this may influence presence beyond the influence of site‐level covariates. The majority (96%) of our surveys were conducted in three years with covariates represented in multiple years and therefore we expect this to have a minimal influence on our results. Nevertheless, an alternative approach to account for such effects would be to incorporate a species‐level linear change or annual random effect in the model.

Conservation implications

We highlight two important conservation implications from this work. First, the concentration of Neotropical migrants at middle elevations in the Colombian Andes overlaps an important agricultural region where there has already been extensive forest loss (Tejedor Garavito 2012). This situation presents challenges for conservation but also opportunities if we can identify management strategies that jointly benefit migratory birds and human livelihoods. A greater emphasis on agricultural practices such as shade‐grown coffee, cacao, and silvopasture offer such opportunities as they have well‐known biodiversity benefits (De Beenhouwer and Aerts 2013, McDermott and Rodewald 2014, Díaz‐Bohórquez et al. 2015) and can improve the quality of agricultural products for landowners (Hernandez‐Aguilera 2018). Second, the diversity of habitats selected by the migrant community suggests there is no one‐size‐fits‐all strategy for conserving migrants overwintering in the Colombian Andes. We argue that those species least tolerant to the more open land covers (e.g., Canada Warbler, Swainson’s Thrush, Acadian Flycatcher) are the ones that need greater attention and, indeed, nearly all are in decline based on long‐term population trends (Sauer 2017, Rosenberg et al. 2019), with wintering ground forest loss being a likely contributor (Wilson et al. 2018). There are limits to forest recovery in these regions given their extensive use as working landscapes. However, a mixed approach combining native forest restoration, where possible, with migrant‐friendly agricultural practices, such as shade coffee, would be a beneficial strategy to help recover these species (González et al. 2020), while also helping to protect resident forest specialists that use these habitats.

Appendix S1

Click here for additional data file.

Data S1

Click here for additional data file.

Data S2

Click here for additional data file.

Data S3

Click here for additional data file.