Mixed-species flocks are a key facilitative interaction for tropical birds. Forest fragmentation leads to species loss and spatial turnover in these flocks, yet it is unknown how these changes to composition influence within-flock species interactions. We used network analysis to characterize flocking interactions along a fragment-size gradient in the Colombian Western Andes. We asked (1) how patch size, edge density, and vegetation structure explained network measures indicative of flock cohesion, (2) whether changes were driven by flocking species turnover or changes to the frequency of species co-occurrence, and (3) whether nuclear species, those that maintain flock stability and cohesion, changed in importance across the gradient. We constructed weighted social networks from flock compositions observed on 500-meter transects, and then calculated global network measures and the centrality of six nuclear species. Patch size and edge density did not correlate with interspecific co-occurrence patterns, but interaction strength increased with canopy height. Flocks contained numerous, weak interactions and there were no flock sub-types, suggesting flock composition was dynamic and unstructured. Several redundant nuclear species were present and varied in importance based on ecological conditions. A chlorospingus (Passerellidae) was most central in old-growth forest, whereas several tanager (Thraupidae) species became more central in smaller fragments and disturbed forest. When partitioning network dissimilarity, we found that 66% of dissimilarity resulted from species turnover, whereas only 34% resulted from changes to species co-occurrence. This finding suggests that coherence of flocking behavior itself is maintained even as extensive species turnover occurs from continuous forest to small fragments.

Study system and sites

We conducted field work in subtropical humid montane forests located in the municipality of El Cairo, Valle del Cauca department in Colombia. The study region is in the Serranía de los Paraguas mountain range, part of the Western Andes, and a center of avian threatened species diversity and endemism in Colombia (Ocampo-Peñuela and Pimm 2014). Andean forests in Colombia are highly fragmented, with only 30% of original forest cover remaining (Etter et al. 2006). Within our focal landscape, we selected eight fragments representing a gradient in patch sizes (range 10 to 170 ha). We stratified forest fragments into large (≥ 100 ha), medium (30-50 ha), and small (≤ 20 ha) size categories and selected a minimum of two replicates of each (Table 1, Jones and Robinson 2020). These sites were in the same altitudinal belt (1900-2200 masl) and matrix type (cattle pasture). Fragments were separated by ≥ 100 meters to minimize cross-patch movement of birds, and all transects in different patches were separated by at least 250 m. We also worked in a non-fragmented reference site (Reserva Natural Comunitária Cerro El Inglés) connected to thousands of hectares of continuous forest.

We only selected fragments with primary forest, though vegetation structure and canopy height varied substantially within and between patches based on varying intensities of selective logging. To capture these within-patch effects of disturbance, we therefore established 500-meter transects through forest interior (N = 14 total transects; Table 1) as our replicate. Where fragment area allowed, we placed one transect in a more disturbed (logged) area and another in a relatively undisturbed site within each fragment. Fragments were privately-owned forests, and we collaborated with a local NGO (Serraniagua; http://www.serraniagua.org) to ensure access. Mixed-species flocks in Colombian sub-montane forests are large and diverse: they contain ~100 participating species and up to 50 individuals per flock (Colorado Zuluaga and Rodewald 2015; Jones and Robinson, 2020). 

Transect surveys for mixed-species flocks

We adapted transect surveys for mixed-species flocks from Goodale et al. (2014), which were conducted within forest fragments from June-August 2017 (boreal migrants absent) and January-March 2018 (boreal migrants present). Both sampling periods correspond to a dry season in the Western Andes, which has a bimodal two dry, two wet seasonality. For each sampling period, we sampled each transect for 2.5 consecutive field days by continuously walking back and forth; we alternated visits to continuous forest, large fragments, and small fragments to avoid a temporal bias in our sampling. Transects were walked slowly by 2-3 observers familiar with all bird species by sight and sound (HHJ present for all surveys), with visits distributed across the morning (7:30-11:30) and evening (15:00-17:30) hours when flocking behavior is most common in Andean forests. When a flock was encountered, we spent up to a maximum of 45 minutes characterizing it with 10x binoculars. We identified flocking species by sight and sound, as flock members often called and sang conspicuously while foraging. At least 5 minutes were spent with each flock, following it if possible. Because detection of bird species in mixed-species flocks is imperfect, we considered the flock composition ‘complete’ when we did not detect new species or individuals over the past 5 minutes of observation; we only used ‘complete’ flock compositions in our analyses. Avian taxonomy of flocking species follows Handbook of Birds of the World Alive (del Hoyo et al. 2020) nomenclature.

Calculation of Landscape-level Variables

            We obtained landscape-level variables for analyses using geographic information software (GIS) analysis in ArcGIS (ArcMap 10.3.1; Esri; Redlands, CA). We quantified landscape composition and configuration by buffering each transect (N = 14) by 1 km. We then calculated measures of landscape composition and configuration using a land cover use categorization made by the Corporación Autónoma Regional del Valle del Cauca, converted to a 25m cell-size raster. To quantify landscape composition, we calculated percentages of forest land use type within each 1 km buffer using the ‘isectpolyrst’ tool in Geospatial Modelling Environment (version 0.7.4.0; Beyer 2015). Following Carrara et al. (2015), we selected percentage of forest cover as a proxy for patch size because some of our transects were located in continuous forest with no patch size measurement. Because the matrix in our landscape consisted of unforested cattle pasture, we feel that this is a good measure of patch area, since there was no other forested habitat in the buffer area; percentage forest was highly correlated with patch area (correlation coefficient = 0.95). We also measured landscape configuration for each transect at the 1-kilometer buffer scale using edge density, or length of all forest edges (in meters) divided by total buffer area (in hectares), as described by Carrara et al. (2015). We did not include other (e.g. 500 m) buffer scales because composition and configuration measures correlated heavily with the 1 km scale.

Vegetation measurements and principal component analysis

            To quantify disturbance to local vegetation at each site, we measured vegetation structure and density along each transect used to sample for flocks. Vegetation measurements were made between June and August 2017 and based on our field observations there appeared to be minimal variation between seasons. We used the methodology of James and Shugart (1970) following the modifications made by Stratford and Stouffer (2013), and further modified it to be used with belt transects. The sampling consisted of two components for each transect: (1) the quantification of canopy cover, ground cover, canopy height, and vertical structure of vegetation using point sampling spaced every 10 meters on the transect and (2) the quantification of shrub, vine, fern, palm, and tree fern and tree density using 3-meter-wide belt transect sampling. Because transects ran along trails, we measured vegetation at least three meters from the trail edge on a randomly selected side for each 100-meter transect segment.

            For the point sampling, we measured eight variables at 10-m intervals, for 50 points per transect. As a measure of vertical vegetation structure along the transect, we noted the presence or absence of live vegetation at five heights: <0.5 m, >0.5–3 m, >3–10 m, >10–20 m, and >20 m. Above 3 meters, we used a rangefinder to determine heights and sighted through a tube with crosshairs while straddling the point. The canopy height at each point was measured using a laser rangefinder (Raider 600 Digital Laser Rangefinder, Redfield Inc. Beaverton, OR) pointed at the highest foliage. Canopy and ground cover were calculated to the nearest 1/8^th of the field of view by sighting through a vertical canopy densiometer (GRS Densiometer, Geographic Resource Solutions, Arcata, CA). For each transect, we averaged values for canopy height, canopy cover, and ground cover, and calculated the proportion of points at which vegetation was present for each height category. For the belt transect sampling, we surveyed vegetation along the same transects and calculated densities for each 100-meter transect interval. We counted all shrubs, vines, ferns, tree ferns, and palms encountered on 1.5 meters to either side. Secondly, we counted all trees (woody vegetation > 2 m in height) within 1.5 meters of the transect and measured their diameter at breast height (dbh). Trees were later categorized into six dbh size classes for analysis: 1-7 cm, 8-15 cm, 16-23 cm, 24-30 cm, 31-50 cm, and > 50 cm. We additionally recorded the largest tree as a measure of degree of logging in each fragment.

            We retained four measures of local vegetation for our analyses. The average canopy height and canopy cover for each transect were directly used for the analyses. Canopy height was strongly correlated with degree of vertical vegetation complexity (i.e. presence of foliage in different height categories), as calculated using the Shannon diversity index on the proportion of points with vegetation present in each of the five height bands for each transect (correlation coefficient = 0.89). We also used principal component analysis (PCA) ordinations of understory vegetation density and density of large-diameter trees, taken from Jones and Robinson (2020), calculated for each transect. We used the first PC axis from each ordination; more negative values of the understory vegetation PC axis indicate higher densities of understory shrubs, vines, palms, ferns, and tree ferns while more positive values of the tree-size PC axis indicate greater densities of large-diameter trees (e.g. 20-50 cm DBH) and therefore reduced selective logging.

Construction, measurement, and analysis of social networks

            We characterized the flocking interactions of the bird community on each 500-meter transect by assembling a social network (Croft et al. 2008). Species move and forage in close association in mixed-species flocks, so we considered two species observed in the same flock to be interacting ecologically (the ‘gambit of the group’ approach; Whitehead and Dufault 1999). We therefore defined each node as an individual species and each edge as a co-occurrence of two species in a flock. All statistical analyses were performed in R (version 3.5.1; R Core Team 2020). We used presence-absence, flock-by-species adjacency matrices derived from field observations to create social networks using the get_network function of the asnipe package (Farine 2013). We constructed one social network for each transect in each sampling period using all flock compositions observed on the transect during that sampling period (range = 7-26 flock compositions per network). However, we did not construct a network for one transect during the boreal winter due to insufficient sample size of flock compositions (N = 27 networks). Because detectability of birds in flocks is high, and associations were likely rarely missed, we used the simple ratio index (SRI), an undirected, weighted measure of association, to calculate an association index for each species pair in the flocking network (Cairns and Schwager 1987).

            For each network, we calculated five global network measures relevant to our questions of interest. We quantified the frequency and strength of species co-occurrences in flocking associations using mean normalized degree, mean weighted degree (hereafter strength), and skewness of strength values for each network. Mean normalized degree represents the average of the number of edges for each node, divided by the total number of nodes minus one. This measure provides an index of network connectedness, the average number of co-occurrence interactions for participating species standardized for network size. Network strength is calculated as the average of all sums of weighted edge lengths for each node in the network, which represents a measure of the consistency of species co-occurrences across flock surveys (compositions). Skewness of strength measures the asymmetry in the distribution of node strength measurements; positive values of skewness indicate that node strength is skewed right (a greater abundance of strength values below the mean), and values approaching zero indicate increasing proportions of high strength nodes (i.e., an approximately unimodal distribution). Mean normalized degree was calculated using the degree function of the igraph package (Csárdi 2019), strength was calculated using the strength function of the igraph package, and skewness was calculated using the skewness function of the moments package (Komsta and Novomestky 2015). Following Mokross et al. (2014), we also used the global clustering coefficient, calculated using the transitivity function of the igraph package, as a measure of flock cohesiveness. Lastly, we quantified the extent to which flocks were divided into sub-types versus homogenous in composition by calculating the modularity of each flocking network. First, we assigned each node (species) to a flocking sub-type using the eigenvector modularity method (Newman 2006) applied through the cluster_leading_eigen function of the igraph package. We then calculated the modularity of this optimal solution using the modularity function of the igraph package.

            In order to understand how nuclear species importance changed across our patch size gradient, we also calculated node-based measures of centrality for six a priori-defined nuclear species. Putative nuclear species were identified based on field observations and descriptions of Andean flocking systems in the literature. We chose Chlorospingus canigularis, Tangara labradorides, T. aurulenta, and Anisognathus somptuosus because they are gregarious species that often join flocks as social groups, a trait often associated with nuclear species (Goodale and Beauchamp 2010), while Myioborus miniatus and Pachysylvia semibrunnea are commonly flocking species that frequently call and sing in the flock, possibly contributing to cohesion. We quantified ‘nuclearity’ as the normalized betweenness centrality of a species in the network, or the number of shortest paths between nodes passing through the node of interest, standardized by the size of the network. We selected betweenness centrality because this measure accounts for not only species with many co-occurrences, but also the ability to connect flock sub-types (network modules), which could be an important “information broker” role (Croft et al. 2008). More broadly, centrality is an appropriate quantification of the importance of the nuclear role because the number and strength of species co-occurrences in flocking networks covary (Srinivasan et al. 2010), and the structural importance of species in flocking social networks is correlated with their functional importance in the flock itself (Sridhar et al. 2013). On one transect, C. canigularis was replaced by C. semifuscus, which plays a similar nuclear role in flocks (Bohórquez 2003) and is ecologically similar (Isler and Isler 1999). We therefore used the centrality value for C. semifuscus on that transect. Centrality was calculated using the betweenness function of the igraph package. 

The following dataset contains three data files: one CSV file containing the presence-absence species composition of each mixed-species flock used to construct 27 social networks for the analysis, one CSV file containing the calculated network measures and associated environmental covariates, and one Word file containing the descriptions of the variables used in the second CSV.