Study area
The study was conducted in the Serranía de los Paraguas region of the Western Andes of southwestern Colombia in the Valle del Cauca Department (Fig. 1). The Serranía encompasses ~ 232,000 ha, composed of hyper-wet Chocóan forest on the west-facing slopes and montane forest more typical of the tropical Andes on the eastern slopes. Our study area was located on the eastern slope in El Cairo municipality (4°45′39″N, 76°13′21″W), between 1250 and 2700 m.a.s.l. This region is part of the Coffee Cultural Landscape of Colombia, a UNESCO World Heritage Site, and recently was categorized by the Colombian National Parks System as a protected area with sustainable use of natural resources (IUCN Category VI; IUCN 2020; UNESCO 2020). This landscape contains a network of ~ 70 privately protected areas (or Natural Reserves of the Civil Society), known as the Reserve Network of Northern Valle del Cauca, of which around half include shade coffee plantations.
The study landscape comprised a mosaic of forest fragments embedded in a matrix of mixed agriculture, primarily shade coffee and cattle pasture, with continuous forest present on the ridgetops. The coffee system supports moderate canopy cover (30–60%) in which more than 80% of the tree species are native, but shade tree diversity is low and mostly planted. Based on photographs acquired by our cameras, some mammals such as tayras (Eira barbara), agoutis (Dasyprocta punctata) and crab-eating foxes (Cerdocyon thous) feed on plantain fruits within coffee plantations. Forested areas in this landscape comprise patches of old secondary forest (typically 20 + years of regrowth) regenerating from pasture or abandoned coffee plantations and, at higher elevations, continuous old-growth forest that connects with the Chocó lowlands to the west (Fig. 1). The topography is steep, with up to 90% slope at higher altitudes.
Focal bird and mammal communities
The potential species pool in this landscape includes approximately 30–35 species of ground-dwelling birds. These bird species, like Andean bird species in general, have “shoestring” ranges within narrow elevational bands, resulting in high beta diversity over the altitudinal gradient (Kattán et al. 2006). Of the expected species at our study elevations, 8 were thought to occur only below 2000 m and 4 only above 2000 m. If there were an altitudinal gradient in species richness in our study landscape, it should therefore be highest at lower elevations in the absence of human disturbance. However, forest fragmentation and disturbance can shift the width of elevational ranges of Andean species (Ocampo-Peñuela and Pimm 2015), and we found some species both above (e.g., Leptotila plumbeiceps and L. verreauxi) and below (e.g., Turdus fuscater) expected elevations at multiple sites. From the species pool for birds, one species (Chestnut wood-quail, Odontophorus hyperythrus) is listed by the IUCN (2020) as Near-Threatened (Appendix A, Table S1), but the ecology and conservation status of many species in our study area are barely known (Greeney et al. 2008).
We expected to detect around 32 mammal species in this altitudinal range of the eastern slope of the Western Andes (Solari et al. 2013). Elevational turnover is less pronounced for mammals than birds in this region (Patterson et al. 1998), and all potential species occur across the entire elevational range of our study area. As with birds, the ecology of many mammals in this region is poorly known (Schipper et al. 2008), but many more species are listed by the IUCN (2020) as of conservation concern: 3 Vulnerable, 5 Near-Threatened, and 2 Data Deficient (Appendix A, Table S2).
Camera trap surveys
We used camera traps to sample birds and mammals in shade coffee plantations and two types of reference forest (forest fragments and continuous old-growth forest) between July 2016 and July 2018. Forest fragments were owned and managed by individual smallholders and generally were part of shade-coffee farms. Therefore, we deployed cameras in fragments and coffee plantations simultaneously at the same sites based on a stratified random sample design by primary land-use type. Each site represented one to four nearby coffee farms containing an agricultural mosaic of shade coffee, forest fragments and pasture (N = 13 sites at 1250–1900 m elevation, black triangles in Fig. 1). Continuous Andean forest sites were located along ridgetops and connected to the continuous Chocoan forest (N = 11 sites at 1600–2700 m, black circles in Fig. 1). These sites were private forested landholdings that were primarily or completely made up of old-growth forests. To ensure that all samples were within the same landscape, all continuous forest sites were located on the east slope of the Andes no more than 8 km from fragment and coffee plantation sites. The minimum distance among all sites was 1 km.
Our sites ranged in area from 3–39 ha for coffee farms and from 20–700 ha for forested landholdings. The number of cameras at each site varied with site area (range = 4–26) for a total of 340 cameras across the 24 sites, and stratified sampling within sites resulted in 57 cameras in shade coffee plantations, 127 in fragmented forest, and 156 in continuous forest. Camera locations ranged in elevation from 1293 to 2518 m (Mean ± SD = 1810 ± 312 m). Cameras were placed a minimum of 100 m apart and deployed for at least one month at each site (Mean ± SD = 39 ± 8 days). We attached cameras to trees ~ 30–40 cm above the ground and did not use baits or lures. All cameras (Browning Trail Cameras, Browning BTC- 5 HD Strike Force) were equipped with motion-activated sensors, and we programmed cameras to take four photographs each time the camera was triggered. We considered consecutive records of the same species to be temporally independent if they were separated in time by at least 24 hours. Cameras within a site were not spatially independent for mammals and some birds because of their large home ranges, and we addressed this lack of independence with a random effect of site in our multi-species occupancy model. Also, we parameterized the detection sub-model to account for variation in camera trapping effort. In our study, occupancy should be interpreted as a measure of habitat use for mammals because species with large home ranges likely moved in and out of our sites during sampling (MacKenzie et al. 2006; Tobler et al. 2015).
For analyses of mammals, we included all species that triggered our cameras except small rodents and marsupials (< 200 grams), which were unidentifiable with cameras. For birds, we only used captures of species that are wholly or partially ground-dwelling, defined a priori based on natural history information on foraging and nesting habits in Birds of the World (Billerman et al. 2020) and expert opinion. Major families included tinamous (Tinamidae), guans (Cracidae), quail (Odontophoridae), doves (Columbidae), antpittas (Grallaridae), and thrushes (Turdidae).
Calculation of landscape variables
To characterize the landscape context of our survey sites, we used remote sensing to classify land use across our study landscape and then calculated landscape-level variables that measured the amount of forest cover (not including shade coffee), degree of isolation, and intensity of human disturbance for each camera location. First, we generated a land use map using a supervised classification with 10-m resolution performed on a Sentinel-3 image acquired through the Google Earth Engine platform (see Appendix B for classification methods). Using this map, we then calculated the following landscape-level variables in ArcMap 10.7 (Esri; Redlands, CA): (1) Percentage of forest, which was the percent of total land cover comprising forest within a 1-km buffer around each camera. We also calculated this variable for a 500-m buffer, but this measure was highly correlated with the 1-km buffer and was not used in analyses. (2) Distance to continuous forest, which we considered a measure of isolation. We calculated the distance to continuous forest as the straight-line distance from each camera location to the nearest edge of a polygon that encompassed continuous old-growth forest. This polygon was built by merging adjacent pixels of the forest class in the land-use map (Appendix B, Figure S1). Strips of forest less than 100 m wide were not included as continuous forest. (3) Human disturbance, represented by an index of human access and influence. We calculated this index using a cost-distance analysis of ease of human access to any grid cell in the study area based on distance to roads, distance to towns (weighted by population), and ease of human movement in the three broad land-use types in our landscape (forest, coffee plantations, and pasture). We used the Summed Point Influence Tool 1.0 beta for these calculations (Fisher and Didier 2012; see Appendix C for detailed methods).
Sampling vegetation and elevation at camera sites
To account for effects of local vegetation on species detectability at camera locations, we sampled vegetation structure around each camera. Variables included: (1) canopy cover calculated as the average of four measures of percent canopy cover taken with a vertical canopy densiometer (Geographic Resources Solution, GRS-Densitometer) at 5 m from each camera in the four cardinal directions, and (2) understory vegetation density calculated by counting the number of 25-cm segments of a 2-m pole obscured by vegetation while standing 10 m away from the pole in front of the camera. We also determined the elevation of each camera with a GPS (Garmin GPSMAP 64s).
Multi-species occupancy model
We used a multi-species occupancy model (Kéry and Royle 2016; Devarajan et al. 2020), implemented in a Bayesian framework, to measure the influence of land-use type (term used for simplicity to represent shade coffee plantations, secondary forest fragments, and older continuous forest) and landscape context on camera-level richness and occupancy of ground-dwelling birds and mammals. In this model, species occupancy is modeled hierarchically to distinguish true absence of a species from non-detection (Royle and Dorazio 2008; Tobler et al. 2015). This framework offers advantages over traditional approaches for inference about species richness and occupancy by accounting for both species-level effects of land-use type and landscape context, as well as aggregated effects of these variables on the full community. This leads to increased precision in estimates of species richness by improving occupancy estimates for all species, including those with low detection rates (Russell et al. 2009).
We constructed the model with the program JAGS (R2Jags: Plummer 2017) implemented in R, which uses Markov Chain Monte Carlo (MCMC) simulation to estimate parameters based on a posterior distribution. We analyzed bird and mammal datasets separately (See model specifications and the final code in Appendix D). We assessed statistical significance of occupancy covariates based on whether 95% credible intervals of beta estimates overlapped 0 and based on visual assessments of violin plots.
Analysis of species richness
We determined species richness for birds and mammals at the camera level and by land-use type for shade coffee and the two types of reference forest. Camera-level richness was calculated within the multi-species occupancy model using a data augmentation approach. We calculated species richness by land-use type using a rarefaction approach because the number of cameras differed across the three land uses. These rarefaction analyses were performed with the package iNEXT (Hsieh et al. 2016), which uses a bootstrap resampling method (bootstrap replicates = 100) and extrapolation sampling curves to estimate an asymptotic species richness. We used the median realized presence/absence (Z) matrix of the multi-species occupancy model, which is an estimate of the observed presence/absence matrix corrected for detection (Dorazio and Royle 2005; Dorazio et al. 2006), as the input for our calculation of richness by land-use type. The Z matrix and the matrix based on naïve occupancy were nearly identical (mismatch between pairs of values = 0.8% for birds, 0.2% for mammals;).
Comparison of species composition between shade coffee and reference forests
To visualize differences in composition of communities among land-use types, we used non-metric multidimensional scaling (NMDS), which measures differences in composition as the variance of distances between the centroids for each land use (i.e., NMDS space). This analysis used a site-by-species presence-absence matrix. To avoid circularity in our analyses, we used the naïve presence-absence matrix and converted it to a pairwise dissimilarity matrix using the Jaccard dissimilarity index in the R package vegan (Oksanen et al. 2016). We subsequently tested for significant differences in community composition across land-use types with a Permutational Multivariate Analysis of Variance (PERMANOVA; Anderson 2017) performed with the function adonis in the vegan package. Because overall differences in composition were significant, we ran a post-hoc permutational pairwise test of significance between land-use types using the RVAideMemoire package (Hervé 2020).