Our research was conducted in Santa Cruz Ocotal Mixe in the state of Oaxaca, Mexico (Fig. 1). Indigenous Mixe communities use the land mainly for coffee production, their principal cash crop. However, the study site is highly modified by human intervention, which has created a heterogeneous landscape with no clear differences among land uses that are interspersed across the territory; this condition promotes sudden habitat changes within short distances. The landscape, therefore, corresponds to a “coffee landscape” sensu Harvey et al. (2021), as it represents a complex mosaic of habitats where coffee is interspersed with other land uses such as secondary forest patches and various types of crops for self-consumption. Shade coffee itself is often diversified with different tree species and crops interspersed between coffee plants, and is the most important crop in this landscape due to its potential economic benefits for farmer communities and dominance in terms of overall area compared to other land uses. Land uses include coffee farms, urban settlements (land dedicated to housing where no coffee or other plantations occur), orchards (bananas, mango, tamarind, guava, avocado trees, or small sugar cane plots), secondary forests (pine-oak, evergreen, and cloud forest remnants), and milpas (plots dedicated to plant corn and beans). Coffee production is dominated by smallholders (0.55–2 ha). Approximately 70% of coffee farms are certified organic. Coffee is grown as a rustic or traditional polyculture shaded system (Moguel and Toledo, 1999) with a canopy dominated mostly by native trees at least 10 m tall and scattered fruit trees. The coffee harvest season extends from December to March, while pruning, weed removal, replanting, and shade regulation is performed at the end of the harvest season in April. Annual crops are produced from June through December across the region.

Fig. 1.

Study landscape. a) Location of the study area highlighted in Oaxaca (dark gray), southern Mexico (16°58'50.8” N, 95°49'25.8” W, 587-1560 masl). b) Location of the 150 point-count stations across the study area.

(2.07MB).

To evaluate the effects of habitat characteristics on the occupancy of insectivorous birds, we randomly established 150 point-count stations with a 30 m fixed radius, separated by a minimum distance of 200 m along established trails in different land management contexts within the study area (i.e., coffee farms, urban settlement, orchards, secondary forests, milpas; Fig. 1) to ensure data independence. In each point-count we measured a total of seven habitat characteristics (i.e., tree cover, shrub cover, herb cover, coffee dominance, epiphytes, tree richness, shrub richness) and one local abiotic factor (elevation) for evaluations of occupancy. Their detailed descriptions are included in Table 1.

Habitat characteristics were selected considering their ecological relevance for regional birds (Ortega-Álvarez et al., 2018) and their importance in coffee farm management (i.e., coffee dominance; Rodrigues et al., 2018). All visual estimations of the habitat characteristics were performed by the same observer to avoid biases. Coffee dominance was particularly important for our study given that coffee plantations were dominant across the landscape, have important implications for the configuration of avian communities, and are directly associated with habitat management (Gregory and Strien, 2010; Philpott and Bichier, 2011; Rodrigues et al., 2018). All habitat characteristics were measured once per season on the first visit to each point-count station; their means and standard deviations are provided in Appendix A in online supplementary material.

Due to the landscape heterogeneity and fuzzy limits among land uses, bird’s occupancy analyses were performed by considering habitat characteristics as explanatory variables instead of land uses. This approach has been previously used by other studies to predict responses of population parameters to habitat gradients (Anderle et al., 2023; Rodrigues et al., 2018), as well as to avoid inconsistent land use classifications (Moguel and Toledo, 1999) and ecologically uninformative results. For our study site, this design further enabled us to: a) identify specific habitat characteristics correlated with bird occupancy, b) assess the effect of habitat characteristics variability on bird occurrence, c) predict occupancy responses to several habitat scenarios across the landscape, including those present in coffee cultivation areas, and d) elaborate more detailed habitat management recommendations.

To sample both resident and winter migratory species across the landscape, two sampling seasons were conducted, one in summer (June and July 2021) and the other in winter (February 2022). Survey months were deliberately selected to avoid regional migratory movements peaks and keep with the closed population assumption associated with single-season occupancy modelling. Sampling was conducted over a 10 min period at each point-count station from dawn until bird activity decreased notably (typically five hours) (Bibby et al., 1998; Ralph et al., 1993). Passing and hovering birds as well as highly mobile species (e.g., raptors, swifts, swallows) were not recorded to avoid double counts when changing point-count stations. Each point-count station was visited twice per season by a group of three different observers trained in the visual and auditory identification of local bird species, which generated a total of 6 observation sessions per station by season (2 visits by 3 observers = 6 observation sessions) to model occupancy. In each visit, the time between each observation (n = 3; i.e., one per observer) in each point-count station varied from 10 min up to 3 days, depending on the proximity of the next station, as suggested by Ruiz-Gutiérrez et al. (2019). This approach enabled us to avoid violations of the occupancy models assumptions of a closed population (Mackenzie, 2006). To prevent biases, the sequence of visits to each point-count station was determined randomly. Single-season occupancy models relied on species detection (detected or not detected) not abundance, and therefore observation sessions made on the same day (n = 6, two per observer per point-count station) did not result in the overestimation of species abundance by double counting of birds, but allowed to establish detection histories required for occupancy estimations (Mackenzie, 2006).

Selected species were subsequently classified into different functional groups of insectivorous birds. We formally constructed such functional groups with a Similarity Profile Analysis using the average cluster method of Euclidean distance from the R package “clustsig” (Whitaker and Christman, 2010). To run this analysis, we considered several species traits defined by field observations and information from the Birds of the World data base (Billerman et al., 2022). Included traits were: main foraging strata (understory or canopy), main foraging strategy (foliage gleaner, bark gleaner, or flycatcher), and body size (tiny <15 g, small 16–35 g, medium 36–60 g, large 61–120 g, and extra-large >120 g). The details of traits classifying species are available in Appendix B.

Prior to constructing the occupancy models, we performed Pearson correlation tests between habitat variables to avoid multicollinearity. When a pair of habitat characteristics were highly correlated (p values above 0.6 or below −0.6), one variable was eliminated keeping the variable which represented fewer challenges for habitat management practices (easier to measure and regulate).

Occupancy modelling was restricted to species from the bird community that could make greater contributions on arthropod pest control across the landscape (Sekercioglu, 2012). We only analyzed species that: a) exhibited a diet of at least 30% insects (Billerman et al., 2022), b) belonged to a taxonomical family previously associated with some type of insect pest control (Michel et al., 2020), and c) were abundant (at least 49 detections) during our surveys. Single-season occupancy models that consider the effect of detectability and habitat covariates on parameter estimation were run for each species using the “unmarked” package of R (Fiske and Chandler, 2011). In particular, we used the “occu” function, which fits the single season occupancy model of MacKenzie et al. (2002) based on zero inflated binomial models (MacKenzie, 2006); logit link was used to model covariates for both occupancy and probability of detection. For a detailed description of occupancy estimation, please see MacKenzie (2006) and MacKenzie et al. (2002). Models were run separately by season, which enabled us to comply with the closed population assumption of the statistical technique and account for seasonal variability in: a) habitat characteristics, b) pests phenology, and c) the functional composition of avian communities. Modeling was performed on species with a minimum of 49 detections in the field to ensure sufficient data for fitting reliable models. To optimize the fit of the models to available data, we estimated probability of detection as a nuisance parameter (Guillera-Arroita et al., 2010), since the performance of the detection method was not of interest for our study. In this way, habitat relationships focused on occupancy estimation while controlling for the effect of imperfect detectability in our models. Occupancy models were constructed for each species considering the additive (no interactions) effect of a maximum of three covariates per model to avoid overparameterization and facilitate the interpretation of results. The three covariates included in the models were chosen by considering all the possible combinations of the total pool of the eight non correlated habitat characteristics measured for this study. A null model where no effect of habitat characteristics is considered was fitted to subsecuently evaluate model fit with different combinations of habitat characteristics.

Model selection was performed using multiple models and an information-theoretic approach. We based model selection on the second order Akaike’s Information Criterion (AICc) because it provides good model selections by adjusting for small sample sizes (Burnham and Anderson, 2002). Finally, we modeled and plotted both positive and negative relations between habitat characteristics and the associated occurrence of each bird species in each functional group (Appendix C). We did not model the response of species richness (total number of species in a functional group) to habitat characteristics because such approach may mask the individual responses of species to habitat traits. To generate response curves of the occupancy probability of each bird species as function of habitat characteristics we used model-averaged estimations.

To synthesize results for different functional groups of insectivorous birds we simultaneously graphed all the curves of bird species from the same functional group responding to a particular habitat characteristic. A minimum of two species per group were required to plot responses to a particular habitat characteristic. Management recommendations were constructed using critical values for each habitat characteristic on our figures, representing the threshold value above which occupancy of most species increased notably. We visually defined critical values as the inflection points where slopes became steeper sensu Ortega-Álvarez et al. (2018), indicating the value in which occupancy probability increased for most species. For cases where there was no inflection point (straighter tendency lines), the average value in the middle of the x axis of the habitat characteristic was considered as the critical value (Ortega-Álvarez et al., 2018). For cases with contrasting positive and negative occupancy probability responses, the critical value was considered as the value of the habitat characteristic where the average of most species’ occupancy was enhanced, as this represents a value where most species still show considerable occurrence. Elevation was not included in management recommendations as it cannot be modified by producers.

ANOVA tests showed that there were no significant differences among land uses for most habitat characteristics, with three exceptions: tree cover, herb cover, and coffee dominance (Appendix D). Tree cover was highest in the secondary forest, followed by coffee farms and orchards in both seasons. Herb cover was highest within the milpa, mainly during winter. Apart from this trend there was minimal variation in the mean values of habitat characteristics across seasons.

We recorded a total of 177 bird species (129 in summer, 145 in winter) belonging to 12 orders and 36 families. Occupancy modelling was performed for 44 insectivorous bird species (Appendix E), of which five were summer migrants, six were winter migrants, and 33 were full year residents for the area. All species were considered habitat generalists as detections of forest specialist species were scarce. The Similarity Profile Analysis identified six functional groups of insectivores at the cutoff of 46% of dissimilarity of Euclidean distance: large foliage gleaners, species between 60 g and 129 g foraging in both understory and canopy; tiny foliage gleaners, species smaller than 15 g foraging on both strata; small and medium foliage gleaners, species between 15 g and 60 g foraging on both strata; small canopy flycatchers, species between 15 g and 35 g foraging on the canopy; bark gleaners, species between 35 g and 120 g foraging in tree bark at both strata; and canopy flycatchers, species between 15 g and 35 g foraging in the canopy (except for Myiopagys viridicata weighing less than 15 g but feeding primarily in the canopy) (Fig. 2).

Fig. 2.

Dendrogram for the classification of insectivorous bird species based on foraging strategies, foraging strata, and body mass. Functional groups are indicated by bird icons and colored boxes at the cutoff of 46% of dissimilitude, matching the color of their functional group name.

(0.65MB).

From the pool of 177 recorded species, we filtered for bird species with the restrictions noted in our data analysis section. After analysis, the occupancy of a total of 34 species in winter, and 22 in summer, was explained by habitat covariates (Appendix F). During summer, shrub cover was the habitat characteristic associated with the occurrence of the highest number of functional groups: tiny foliage gleaners, small and medium foliage gleaners, small canopy flycatchers, and bark gleaners (62.5% of species with a positive response and 37.5% of species with a negative response) (Table 2).

As shown in Fig. 3, occupancy of a particular number of species from each functional group was influenced by habitat characteristics. For instance, shrub cover was positively associated with all tiny, small, and medium foliage gleaner species, with the only exception being Troglodytes aedon, a tiny foliage gleaner adapted to rocky habitats with little to no vegetation. Tree cover was negatively associated with 7 out of 9 species from all three groups of foliage gleaners that responded to this trait, both species of large foliage gleaners, three small and medium foliage gleaners and two tiny foliage gleaners (Fig. 3). Herb cover was correlated both positively and negatively with the occupancy of large foliage gleaners (1 out of 2 species, respectively) and tiny foliage gleaners (1 out of 2 species). Coffee dominance was positively associated with 2 out of 3 tiny foliage gleaners, excluding Pheugopedius maculipectus.

Fig. 3.

Projected response of the occupancy probability of each insectivorous bird functional group to habitat characteristics during summer and winter. Each line represents a different insectivorous bird species. The line type, the color of the line, and the bird icon depicts the functional group of the species.

(0.55MB).

During the winter, tree cover negatively affected small and medium foliage gleaners (3 out of 3 affected species), and tiny foliage gleaners (1 out of 2 affected species) with the exception being Basileuterus culicivorus (Fig. 3). Herb cover was positively associated with small and medium (3 out of 3 affected species), and large foliage gleaners (3 out of 3 affected species). Coffee dominance had a positive relationship with three functional groups: bark gleaners (4 out of 4 affected species), small canopy flycatchers (4 out of 4 affected species), and tiny foliage gleaners (1 out of 2 affected species), with the exception of Basileuterus rufifrons. Epiphytes were positively correlated with the occurrence of every species affected by this trait from a total of five functional groups: tiny foliage gleaners (3 species), small canopy flycatchers (2 species), large foliage gleaners (3 species), canopy flycatchers (3 species), and bark gleaners (1 species) (Appendix E). Finally, tree richness positively affected the occupancy of large foliage gleaners (2 out of 2 species), whereas shrub richness was positively correlated with the occupancy of small and medium foliage gleaners (2 out of 2 affected species).

The estimated critical values where most functional groups increased their occupancy probability for each habitat characteristic are shown in Appendix A. Tree cover only exceeded its proposed critical value within coffee farms and secondary forests in both seasons. Shrub cover was higher than its respective critical value for both seasons in milpas and orchards, but only higher in coffee farms during summer. Herb cover only exceeded its proposed critical value within milpas during winter. Epiphytes never reached or exceeded their proposed critical values; highest values were found in orchards during summer and in secondary forest during winter.

An increased number and diversity of insectivorous birds can enhance pest control services in coffee farms, specially by foliage gleaning species from the Parulidae family and flycatching birds from the Tyrannidae family (Kellermann et al., 2008; Martínez-Salinas et al., 2016; Sherry et al., 2016). Given that occupancy should be positively correlated with abundance (MacKenzie, 2006), management recommendations based on the thresholds of specific habitat characteristics may be useful for simultaneously increasing insectivores’ occupancy and pest control service across the study landscape. Nevertheless, we recognize that pest control is a complex ecosystem service that can involve multi-trophic interactions between service providers and agricultural pests that might not always result in higher yields by enhancing habitat conditions for birds in the landscape (Díaz-Siefer et al., 2022). For example, depending on the study system, a trade-off between pest predation and pest predation alleviation due to birds consuming other non-pest, pollinating, or even other pest controlling invertebrates, could result in a dampened ecosystem service (Garcia et al., 2021). In contrast, for coffee agroecosystems, a complementary effect has been observed between bird pest control and pollination services, as both have shown to increase the quantity and quality of coffee fruits respectively without negatively affecting one another (Classen et al., 2014) and even creating positive synergies (Martínez-Salinas et al., 2022). Further research on this topic in the same region may help shed light on multi-trophic interactions and their effect on pest control service.

Given that coffee farming in Latin America is usually practiced in diverse and heterogeneous landscapes (Harvey et al., 2021), our findings suggest that the use of habitat characteristics can be a practical approach to determine avian-habitat relationships. As birds do not exclusively occupy coffee farms, it is important to measure these characteristics throughout the landscape to identify their associated values and possible levels of specific ecosystem service provisioning (e.g., pest control). In this sense, we identified specific critical values of habitat characteristics that might contribute to increase the occupancy probabilities of different insectivorous birds and ultimately enhance their potential pest control service. Specific management recommendations related with such critical values include maintaining a tree cover around 45% inside coffee farms by favoring tree pruning over removal, and planting more native tree species or allowing them to establish naturally. We realize that the relation between shade and coffee yields is important for producers, as it can vary between coffee cultivars (Koutouleas et al., 2022). However, we expect that our proposed level of 45% can ultimately be beneficial for coffee yield production, as this tree cover can favor pest control provision by increasing the occupancy probabilities of tiny, and small and medium foliage gleaning insectivores. Based on our results, we recommend promoting or propagating epiphytes’ cover above 27% on shade coffee trees, suggesting an important threshold versus previous studies by Cruz-Angón and Greenberg (2005) contrasting epiphytes presence-absence only. We also recommend preserving a native shrub cover of 48%, and herb cover of 50% within farms or along their borders, as such values may be particularly important in benefiting foliage gleaners and all flycatcher functional groups. Avoiding epiphyte removal, identifying the epiphyte and shrub species native to the area, and cultivating them in or around coffee farms could be important in this regard. Keeping a diverse tree and shrub richness of 13 and 17 species respectively might also be important for small, medium, and large foliage gleaning insectivores. The presence of these insectivores could be promoted by diversifying crops inside farms or establishing living fences. Although marked dominance by coffee may favor the occupancy of tiny foliage gleaners, small canopy flycatchers, and bark gleaners, dramatic increases should be avoided as they typically require reductions in the number or diversity of shade trees that are also important factors in insectivore abundance and are typically associated with an increased use of agrochemicals.

Our results regarding the medium levels of tree cover that increase the occupancy probability of most foliage gleaners, support the current levels of tree cover (e.g., 40%) required by most conservation-based certifications (Rainforest Alliance, Smithsonian Bird Friendly). However, shrub, herb, and epiphytes cover, as well as tree and shrub species richness should be considered when seeking to maximize ecosystem functionality. The absence of unmanaged forests in the landscape should also be a major concern. Logging in adjacent plots of forest may lead to a reduction of insectivorous birds abundance (Vergara et al., 2021) and to the loss of forest-dependent species, which in turn could ultimately increase insect herbivory in a trophic cascade (Peter et al., 2015).