Habitat alteration modifies the structure and function of mixed‐species flocks in an Andean landscape

Understanding the responses of species interactions to land‐use change is a key challenge in ecology and conservation. In the tropical Andes, a large proportion of birds interact with other bird species and create networks of mixed‐species flocks; however, the widespread effects of land‐use change in the region raise concerns about the alteration of these interactions. Here we explored how structural (i.e., species richness and size) and functional (i.e., network connectivity indices) parameters of mixed flocks change along a gradient from open vegetation to native forest habitats. This study was conducted in the southern Andes of Ecuador, where we sampled mixed‐species flocks in 12 areas for five consecutive months, and applied network analysis tools to characterize structural parameters (observed richness, flock size, flock encounter rate, composition) and functional parameters (degree, weighted degree, cluster, and skewness) of flocks. We found that the network indices, degree and weighted degree of species, were positively associated with areas of dense vegetation and native trees. These two network indices indicated that mixed flocks occupying native forest have higher connectivity, a factor that can promote the stability of flocks. The structural parameters of flocks, species richness, and abundance, were not associated with the vegetation structure. Overall, our results show that functional parameters of flocks can be more sensitive than structural parameters to the effects of land‐use change in our study area. We demonstrate the importance of combining the analysis of structure and functionality of flocks to reveal the effects of land‐use change. Lastly, our results highlight the importance of mature native forest for the conservation of species and the maintenance of non‐trophic interactions.


| INTRODUC TI ON
Land-use change affects the abundance of populations, distribution of species, and recent findings indicate that it can also impact species interactions (Alroy, 2017;Barnes et al., 2017;Jones & Robinson, 2020;Karp et al., 2012). Although it is known that trophic interactions like pollination and seed dispersal can be negatively affected by land-use change (Campbell et al., 2018;Maruyama et al., 2019;Quitián et al., 2018), resulting in extinction cascades (Montoya et al., 2006;Morris, 2010;Säterberg et al., 2013), the effects on non-trophic interactions, such as mixed-species flocks of birds, are often overlooked. However, these associations can influence community structure and ecosystem functioning (Goodale et al., 2010;Traill et al., 2010). Approximately 38% of bird species in tropical ecosystems are known to participate in mixed-species flocks (Zou et al., 2018), and most species within flocks participate in seed dispersal and control of insect populations, two key processes for ecosystem functioning (Munoz, 2016;Srinivasan et al., 2012;Wunderle, 1997). The recognition of the importance of mixed-species flocks has given rise to concerns about the breakdown of these interspecific interactions due to anthropogenic disturbance in tropical areas (Zou et al., 2018).
Mixed-species flocks of birds can be particularly sensitive to anthropogenic habitat disturbance (Borah et al., 2018;Goodale et al., 2015;Valiente-Banuet et al., 2015). Ultimately, bird species join flocks due to positive fitness effects related to complementary vigilance of predators and/or increased foraging efficiency (Sridhar et al., 2009(Sridhar et al., , 2012Terborgh, 1990); however, land-use change alters habitat resources and environmental conditions, which may result in an overall change in the benefits of flocking and the propensity to join flocks (Mammides et al., 2015;Zhang et al., 2013). As a result, common responses of mixed-species flocks to habitat disturbance are a decrease in the species richness and size of flocks (Knowlton & Graham, 2011;Maldonado-Coelho & Marini, 2000;Zuluaga & Rodewald, 2015).
Traditionally, studies on mixed-species flocks have focused on describing the response of flock structural parameters (i.e., richness, abundance, and composition) to environmental gradients; however, species within flocks are part of functional interactions that can be better studied using network analysis tools (Farine et al., 2012;Mammides et al., 2018;Mokross et al., 2014;Montaño-Centellas, 2020;Wey et al., 2008). In mixed-species flock networks, species are treated as nodes, and the connections among species are links based on species co-occurrences (Farine & Whitehead, 2015;Kay et al., 2018). Network analysis provides measures of direct and indirect interactions in flocks at the level of individuals, species, and flocks, which can be linked to emergent properties of ecological functioning such as stability (Farine & Whitehead, 2015). For instance, the exchange of information is a currency that results in increased fitness gains among members of mixed-species flocks (Sridhar et al., 2009); thus, any disruption in the flow of information among flock members will ultimately alter flock stability. Network metrics give insight into the flow of information in networks by measuring connectivity levels among interacting members and thus characterize the stability of mixed-species flocks and their responses to environmental change (Goodale et al., 2010;Mokross et al., 2014;Montaño-Centellas, 2020).
Network analysis tools have revealed a decrease in the frequency of interactions among species and a reduced cohesion of flocks in disturbed areas (Mammides et al., 2018;Mokross et al., 2014).
Apparently, species that join flocks in disturbed environments do so more opportunistically, resulting in less connected and stable flocks.
However, the use of an interaction networks framework to study responses of mixed flocks to disturbance is not new, but there is still a need for more knowledge from tropical areas. For example, there are not yet studies from the tropical Andes, an area that harbors one of the largest concentrations of species of birds in the world (Reid, 1998), many of which are facing conservation problems driven primarily by land-use change (Brooks et al., 2002;Myers et al., 2000).

| Study area
This study was conducted in the southern Andes of Ecuador, in Cajas National Park and surrounding areas of Azuay province. This region is a global priority area for the conservation of biodiversity (Astudillo et al., 2015;BirdLife International, 2014;Stattersfield et al., 2012).
Cajas has been affected by anthropogenic activities for centuries, and the current landscape is composed of a mosaic of vegetation types in different successional stages (MAE, 2011).
We worked in three different valleys: Mazán, Llaviuco, and Sayausí ( Figure 1). Mazán is a U-shaped valley, with a strict protection status that only allows research activities. Mazán is covered by primary montane forest and secondary forest, while the lower areas of the valley are dominated by introduced tree species, characterized by mixed stands of mature Eucalyptus sp., Pinus radiata, and native species in the understory. Llaviuco is also a U-shaped valley, located within Cajas National Park. This is a recreational area, with activities such as hiking and fishing restricted to select locations.
Most of the forest in Llaviuco was cleared for cattle ranching, but, since being incorporated into Cajas National Park in 1996 and the removal of the cattle, native and non-native vegetation has started to naturally reestablish. Sayausí valley is an agricultural valley that has lost much of its natural vegetation and is currently dominated mainly by pastures, with forest patches of Eucalyptus trees, native vegetation remnants in inaccessible areas, and live fences of native species serving as boundaries between livestock areas.

| Study design and sampling
In each of the three valleys, we delineated four transects of 500 m, each separated by a distance of 250 m in order to avoid spatial overlap among flocks (Zuluaga & Rodewald, 2015). Due to the heterogeneity of the vegetation in the study area, instead of using categorical classifications of vegetation types, we decided to obtain quantitative measures of the vegetation structure. Along each transect, we systematically placed six 5 m radius circular plots every 100 m; data from the plots were pooled to quantify five vegetation characteristics that are known to influence bird diversity (Schulz, 2009;Zurita et al., 2006): number of shrubs, number of small and large trees, number of exotic trees, canopy cover, and canopy height. In each plot, we counted the number of shrubs (woody plants <3 cm in diameter at breast height (DBH) that touched the extended arms of a person walking along the two diagonals located in the center of plot directed north, south, east, and west transect. In the same transect, we estimated the canopy cover at 1 m intervals. Canopy cover was quantified by looking directly up at the canopy through a 10 cm diameter tube, which had been divided in four equal parts; by counting the number of quarters covered by vegetation, we estimated canopy cover on a scale from 0 to 4, where 0 indicated no coverage and 4 maximum coverage. The number of trees inside the plot was counted and classified as small (DBH: 3 to 23 cm) or large (DBH >24 cm). Canopy height was evaluated by estimating the height of the 20 highest trees in the plot. Finally, we counted the number of exotic trees, mainly Eucalyptus globulus, present in each plot. All these measurements were summarized for each transect by calculating the means of the six circular plots per transect.
Each transect was sampled for flocks once in the morning, between 6:00-10:00 and once in the afternoon from 15:00-17:00, during four consecutive days, which corresponds to one sampling period. The starting time of sampling in each transect varied within a sampling period with the intention to capture the variation in the daily activity of birds across transects. During sampling, two observers walked each transect for 30 to 40 min (depending on flock detection), making short stops of three to five mins every 100 m. We used a fixed width of 30 m on each side of the transect to search for flocks (Rosenstock et al., 2002;Verner, 1985), where a flock was defined as individuals of two or more species observed traveling and foraging together for more than 5 mins (Latta & Wunderle, 1996). When a flock was detected, it was observed and followed, even outside the fixed width of the transect, for a maximum of 15 min. For each flock detected, we recorded the identity and number of individuals of each species present. In total, we carried out five sampling periods on each transect from June to October 2018.

| Statistical analyses
To summarize the vegetation characteristics of each transect, we used a principal component analysis (PCA). We employed this method because we wanted a condensed description of the vegetation based on all five measured variables; moreover, since flocks are formed by species with multiple ecological requirements (Morse, 1970;Powell, 1985), a combined descriptor of the vegetation is a proper way for establishing ecological responses at the flock level.
The PCA was based on a correlation matrix to consider that original variables have different units (Table S1). The PCA generates a new set of orthogonal variables composed of linear combinations of the original variables, which contain the greatest variance among the original variables (Jolliffe & Cadima, 2016). To perform this analysis, we used the princomp function in the R environment (Kassambara, 2017;R Core Team, 2018).
Based on a screen plot of the broken-stick method (Jackson, 1993), we retained the first three PCA components to characterize the structure of the vegetation in each transect. These PCA components represented 86% of accumulated variance of all the measured  Figure   S1). Hence, we consider PCI to describe a gradient of native shrubs and small trees to an increasing presence of exotic trees; PC II to represent a gradient of short vegetation and shrubs to large native trees and high canopy cover, and PCIII to describe a gradient of open vegetation to closed vegetation with increasingly tall trees.
We obtained structural (flock size, richness, encounter rate, and composition) and functional (network metrics) parameters of flocks for each sampling period on each transect. Flock size (number of individuals in a flock) and observed richness (number of species in a flock) and encounter rate (number of flocks per transect and sampling period) were obtained for each transect per sampling period. Weighted interaction networks were built for each transect during each sampling period (Farine & Whitehead, 2015;Wey et al., 2008). The interaction networks were constructed using the "group gambit" method (Franks et al., 2010), which defines all the individuals of a group of animals observed at any given time as associated (Whitehead & Dufault, 1999). This method assumes interactions among all the species that were part of a flock and is a common method employed in networks of mixed-species flocks (Mokross et al., 2014;Montaño-Centellas, 2020;Zhou et al., 2019). To build weighted networks, we used the frequency of co-occurrence of species in flocks within a transect during each sampling period. We obtained five weighted and adjacent symmetric interaction networks for each transect, for an overall of 59 networks (one transect could not be visited during a sampling period for logistical reasons) in this study.
We obtained network metrics that provide important information about the functional connectivity of flocks at the level of species and flocks (Farine & Whitehead, 2015;Wey et al., 2008), described in Table 2: degree, weighted degree, cluster coefficient, and skewness. To determine if observed networks could be differentiated from random associations, we created null models using a node permutation algorithm (Network permutations, Figure S3) (Farine, 2013(Farine, , 2017a). This null model redistributes all the attributes of the nodes in the network, while maintaining the same number of nodes (Farine, 2017a;Farine & Whitehead, 2015). To calculate the node metrics, we ran 1000 permutations using the R packages Igraph (Csardi & Nepusz, 2006) and Asnipe (Farine, 2017b), and visualized the networks with the package Network (Butts et al., 2019).
The number of animal species detected in field studies can vary across habitats with different vegetation structures (Boulinier et al., 1998;Gu & Swihart, 2004); it could be easier to detect birds on open habitats compared to habitats with a denser vegetation. Therefore, we evaluated if sampling completeness in species richness was influenced by the vegetation structure. A sampling completeness index was calculated as the percentage of species observed in a flock over the predicted species richness obtain by the Chao I index (Chao, 1984;Gotelli & Colwell, 2001;Kéry & Schmid, 2005). Then, we constructed linear models using sampling completeness as response variable and the PCAs of vegetation structure as predictor variables.

Degree Species
The sum of how many species interact with a particular species.

Higher degree indicates more interspecific interactions
Weighted degree Species It represents the association rate per species. The number of times (frequency) a species had an interaction with another particular species.
Higher weighted degree indicates a higher frequency of interactions.

Cluster coefficient Network
The proportion of species connected with the focal species and that are interconnected.
Higher cluster coefficient indicates a greater grouping of species in the network.

Skewness Network
Quantifies the level of asymmetry in the frequency distribution of the number of connections between species within a network.
Indicates that species have similar or different degrees.
Note: More details about them can be found in Farine and Whitehead (2015) and Wey et al. (2008).
To explore how the structure of the vegetation in the sampled area influenced the structure and function of the flocks, we used linear mixed models and generalized linear mixed models (Harrison et al., 2018). The vegetation structure represented by the PCAs was used as predictor variables, and the different variables of flock structure (observed richness, size, and encounter rate) and function (degree, weighted degree, cluster coefficient, and skewness) as response variables in univariate models. For abundance, richness, degree, and weighted degree, we used the function glmer with a Poisson distribution and log-linear link (McCullagh & Nelder, 1989).
For encounter rate, cluster coefficient, and skewness, we used the function lmer (Bates et al., 2014;Harrison et al., 2018); the functions glmer and lmer are part of the Lme4 package (Bates et al., 2014).
Given the hierarchical spatial structure of our study design (transects located within valleys), we nested transects within valleys as a random factor in all the models (Schielzeth & Nakagawa, 2013).
We also used the sampling period as a random factor to account for potential temporal autocorrelation among samples. Statistical significance in the models was based on p-values (<0.05) and coefficients with 95% confidence intervals that did not overlap with zero.
A diagnostic assessment of the models was performed by estimating overdispersion (Table S4) and by visually inspecting the distribution of residuals ( Figure S3).
To test for the effects of vegetation structure on flock composition, we used a multivariate generalized linear model (GLM) implemented in the package mvabund . This is a test fits a GLM for each species and accounts for a better mean-variance relationship than distance-based ordination methods . It allows for testing responses at the community and species level. In these models, we used PCI, PCII, and PCII as predictor variables, fitting a negative binomial distribution for each species.
A Wald test statistic and p-values were used to test for communitylevel effects of the predictor variables, and adjusted p-values were calculated for each species recorded in the flocks. All analyses were performed in the R environment (R Core Team, 2018).

| RE SULTS
We  (Table 3; Table S2). Sampling completeness in species richness was not associated with any of the variables that describe the vegetation structure (PCI, PCII, PCII); thus, the number of species detected in flocks across samples was likely not influenced by differences in the vegetation structure (Table S5).
There were different associations between the structural flock metrics and vegetation structure (Figure 2; Table S1). The observed richness and size were not associated with vegetation structure, but the encounter rate decreased in areas with tall exotic trees (PCI) (Figure 2). Flock species composition was significantly influenced by the abundance of exotic trees (PCI) and native trees (PCII) ( Table 4). Two species responded significantly to an increasing presence of exotic trees: Cranioleuca antisiensis responded positively, while Anisognathus igniventris responded negatively (Table 4). Likewise, two species responded to an increasing presence of large native trees with high canopy cover: Myiothlypis coronate responded positively and Anairetes parulus negatively (Table 4)

| DISCUSS ION
This study shows that mixed-species flock can respond to changes in vegetation structure associated with land-use change in the tropical

TA B L E 3
The frequency of flock participation (how many times a species was observed in a flock) of the 10 most abundant species observed in the mixed-species flocks of birds in the tropical Andes of Azuay province, Ecuador Andes. However, the only structural parameter of flocks that was as-   (Thiollay, 1999).
Therefore, the continued participation of some species in flocking in open areas may result in flocks with tighter interspecific associations (Borah et al., 2018).
The encounter rate of flocks was the only structural parameter of flocks that was associated with the structure of the vegetation, showing a decreased in areas with characterized by the presence of the exotic Eucalyptus trees. The encounter rate of flocks is a parameter that can be sensitive to anthropogenic disturbance . Eucalyptus trees are widespread in the Andes (Aide et al., 2019;Farley, 2007), and they could have negative effects on the diversity of plants and animals (García-Suabita et al., 2019;Gardner et al., 2007;Latta et al., 2011;Rios et al., 2015). Here we show that the presence of Eucalyptus trees might not be favorable to the formation of mixed-species flocks.
We also found that the identity of species participating in flocks can be influenced by vegetation structure and the presence of exotic trees. Changes in species composition of flocks related to land-use change have been described in multiple ecosystems (Brandt et al., 2009;Jones & Robinson, 2020;Montaño-Centellas, 2020;Thomson & Ferguson, 2007). In our study, the variation in species composition  ) could be sensitive to habitat disturbance (Mammides et al., 2015), and the decline or loss of those species can reduce flocking propensity and even cause flock breakdown (Mammides et al., 2018;Rutt et al., 2020;Sridhar & Sankar, 2008).
Disentangling which of these mechanisms is producing the observed changes in flock characteristics is complex, but obtaining comple-    Table S3) A potential limitation of our study is that we did not consider differences in the detectability of individuals or species within flocks in sampling areas with different vegetation structure. That problem could potentially produce biased estimates of some of the structural and functional parameters of flocks.
However, sampling completeness was similar across areas (Table   S5). Our sampling protocol, where we had two observers simultaneously sampling flocks within a narrow width fixed transect, likely decreased the influence of the variation in detectability of species across transects (Rosenstock et al., 2002;Verner, 1985).
Moreover, statistical simulations of networks have found that functional parameters measured by network analysis can be highly robust to missing individuals (Farine & Whitehead, 2015).
Therefore, the variation in structural and functional parameters observed should be mainly driven by the biological responses of species and flocks.
The formation and maintenance of mixed-species flocks are key for the conservation of bird diversity in the tropical Andes (Jankowski et al., 2010;Rutt et al., 2020), and our study shows that mixed-species flocks can be sensitive to the effects of anthropogenic disturbance. We highlight the value of using network metric tools to explore the function of mixed-species flocks and their responses to disturbance. By complementing the information provided with structural (e.g., number of species and individuals) and functional (e.g., number and type of interactions) parameters, it is possible to have a thorough picture of the responses of mixed-species flocks to the current impacts of global change.

ACK N OWLED G EM ENTS
BT and BV were funded by Vicerrectorado de investigaciones-Universidad del Azuay, grant number: 2019-0096 and also the Secretaría de Educación Superior, Ciencia, Tecnología e Innovación (SENESCYT). We also thank UCUENCA EP for the financial administration of the project. We also acknowledge Michelle Armijos and the intern from Universidad del Azuay who helped to collect data in the field.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in the Dryad Digital Repository: https://doi.org/10.5061/dryad. xsj3t x9dx (Vásquez-Ávila et al., 2021).