Hummingbird–Plant Interactions Are More Specialized in Forest Compared to Co ﬀ ee Plantations

: Deforestation transforms habitats, displacing vertebrates and the other dimensions of biodiversity they support through their interactions. Few empirical studies have quantiﬁed the e ﬀ ect deforestation has on vertebrate–pollinator interaction networks. Here we quantify how hummingbird–plant networks change in relation to hummingbird diversity across a deforestation gradient. We found that, overall, hummingbird–plant interactions were signiﬁcantly more specialized in forests and specialized interactions decayed rapidly with the loss of tree cover at small spatial scales. Hummingbird species interaction specialization was also higher in forest habitats compared to co ﬀ ee plantations, but we found no support for a morphological hummingbird trait that predicted interaction specialization or forest dependence. Finally, we developed spatially explicit models for quantifying impacts of land-use decisions on hummingbird species and the biodiversity they support. These tools can be used to identify and prioritize important habitats for conservation activities, like creating new protected areas and improving agricultural lands for biodiversity.


Introduction
Conversion of natural habitat into land used for agricultural production is a major threat to biodiversity, particularly in tropical forests [1][2][3]. Human-made habitats of cropland and pasture make up~40% of Earth's ice-free terrestrial area, and they tend to occupy the most fertile and biodiverse parts of a landscape [4,5]. Deforestation is typically driven by economic opportunities motivating landowners to make way for agriculture, and rates of deforestation continue to climb [3,6]. Despite the reality that agricultural habitats are ubiquitous in the Anthropocene, their capacity to support biodiversity is rarely explored, but studies suggest that some agricultural landscapes have conservation potential. For example, a recent study reported that agricultural landscapes across Latin America support substantial biodiversity, as 45% of 732 vertebrate species were recorded in agricultural plots, 100 of which were found exclusively in agriculture [7]. Given the limitations of expanding protected areas in many regions, a deeper exploration of the biodiversity in agricultural habitats is warranted, especially in light of recent recommendations to intensify agriculture under the optimistic assumption that it will coincide with the sparing land of equal size, but higher quality [8].
A major knowledge gap in quantifying biodiversity in agricultural lands is the dearth of studies on species interactions at the community scale, especially those that compare species interactions in agricultural habitats to those found in remnant forest habitats or reserves. Species interaction networks quantify the species interactions within an ecological community and provide insight into community assembly, biodiversity persistence, ecosystem functioning, and species' resource use [9][10][11].
Of particular interest in the Neotropics are hummingbird-plant interactions, as 10-30% of plant species in Neotropical forests are primarily hummingbird pollinated [12][13][14][15]. Changes in hummingbird-plant interactions can therefore alter plant population dynamics with long-term consequences for forest composition [16,17]. While we know that deforestation can impact hummingbird diversity [14,18,19], given the challenge of collecting species interaction data in diverse tropical regions, empirical data on deforestation's effects on hummingbird-plant interactions remain scarce [20]. Furthermore, the capacity for agricultural lands to support hummingbird diversity is unknown, and so are the consequences of intensifying agriculture by removing remnant forests and clusters of trees in the tropical regions that grow the world's coffee, chocolate, and other products.
Species interaction networks can also reveal how communities change and which species are winners and losers when native forest is converted to agriculture. For instance, specialist species (those with few interaction partners or highly selective interaction partner choice) are typically more vulnerable to environmental change, and thus less common in human-made habitats [9,21]. Greater interaction partner specificity indicates greater niche partitioning within a community, meaning that specialized species may reduce interspecific competition, allowing for greater biodiversity maintenance [22][23][24]. Hummingbird morphological traits, such as bill length and curvature, are hypothesized to play a role in hummingbird interaction specialization and, consequently, hummingbird response to development of agricultural habitats [19,20,25,26]. Despite the vulnerability of specialized interactions, many insect pollinator species have been found to be quite labile in their interaction partner selection [10,[27][28][29]. Moreover, [20] found that the relationship between hummingbird morphological traits and interaction specialization was dependent on resource availability. Empirical evidence of hummingbird-plant interactions inclusive of both agricultural and forested habitats will help resolve the link between morphology, interaction specialization, and vulnerability to deforestation and agricultural intensification.
Here we used a dataset of 3614 hummingbird-plant interactions collected along a deforestation gradient in southern Costa Rica to explore changes in interaction behavior and network structure. The 19 study sites were located in forest reserves, forest fragments, and coffee plantations, but deforestation was also measured at a fine-scale by mapping the locations of individual trees throughout the landscape, allowing us to quantify the contributions of small-scale diversifying features in agricultural habitats thought to promote hummingbird diversity and support hummingbird-plant interactions. In addition, using field and museum specimens, we measured four hummingbird traits thought to be associated with hummingbird specialization. Overall, we aimed to answer the following questions: (i) How does conversion of native forest into agricultural habitat influence hummingbird diversity, community composition, and hummingbird-plant network structure? And (ii) what is the relationship between hummingbird interaction specialization and morphological traits, and do they predict hummingbird species dependence on forest habitat or confer an advantage for colonizing coffee plantations?

Methods
We sampled hummingbird biodiversity and quantified tree cover across a coffee producing region in Costa Rica (8 • 47 N, 82 • 57 W, 800-1400 m elevation, Table S8). Study sites were situated along a gradient of deforestation that can be lumped into three major categories: forest reserves, forest fragments, and coffee plantations. Forest fragments varied by size and quality but consisted of stands of trees with closed canopies overhead (mean ± SD = 86% ± 10% permanent canopy cover). Coffee plantations in the region do not adhere to any convention of shade-grown or organic but are diversified in structure by having a mean ± SD of 28% ± 26% seasonal canopy cover mostly consisting of Erythrina peoppigiana trees and Musa species (banana and plantain plants).
The region was entirely forested prior to 1947, but today it is characterized by a mosaic of pastures, croplands, coffee plantations, forest fragments, and forest reserves of various sizes and qualities that were probably created during peak deforestation between 1960 and 1980 [30]. Today, exactly 40% the 654 km 2 of the land outside the protected areas in the agricultural ecosystem is covered by trees [7]. Tree cover was quantified across the entire region using Google Earth to manually digitize all trees on the landscape-and attempted to include all scattered trees and stands regardless of size or quality [31].
Hummingbird sampling occurred over 5 years of mist netting at 19 sites in a Costa Rican agricultural ecosystem, including 4 sites located inside forest reserves, 9 sites in forest fragments, and 6 sites in coffee plantations ( Figure S1). Hummingbird sampling included a forest reserve site inside La Amistad Biosphere Reserve (4010 km 2 ) located~25 km away from the Las Cruces Reserve. Manipulations of hummingbirds were conducted and approved under Stanford University's Institutional Animal Care and Use Committee (IACUC), Assurance Number: A3213-01, Protocol: 26920 approved in 2006. We supplemented our hummingbird species richness and abundance data using published results from hummingbirds sampled from 14 different forest fragment sites operated nearby [19]. In [19] hummingbirds were sampled using similar protocols to us over the same time period and in the same study area (average distance between sites, excluding one site located in the La Amistad Biosphere Reserve, ± SD = 1.9 ± 0.95 km; Figure S1). Therefore, we model hummingbird species richness and abundance using data from 4 sites located inside forest reserves, 23 sites in forest fragments (9 from our sampling and 14 from data published by [19]), and 6 sites in coffee plantations. Data from [19] were only used in analysis of hummingbird species richness and abundance in the agricultural ecosystem.
We measured hummingbird diversity in the agricultural ecosystem by modeling hummingbird species richness using Chao species richness estimates that account for detection bias [32,33], including hummingbird species richness for 14 forest fragment sites reported by [19]. Hummingbird abundance based on unique hummingbird captures and abundance values from [19] was doubled to account for 50% less sampling based on their use of 10 versus 20 mist nets to capture birds. Hummingbird capture data from [19] were not used in analyses of community composition, which was calculated using Bray-Curtis dissimilarity coefficients and tested for differences among site categories using a permutational MANOVA with the function adonis in the R package vegan version 2.4.4 [34].
Tree cover was used to conduct a multiscale analysis to select a spatial scale of effect [35] that best explained hummingbird species richness based on the proportion of tree cover inside a specified spatial scale. For the multiscale analysis, models were constructed using a single value for the proportion of tree cover ranging from 10 to 1000 m from mist net locations at each site. Models that used site identity as a discrete category of habitat type (i.e., forest reserves, forest fragments, and coffee plantations) and binary categories of habitat type (i.e., forest sites versus coffee plantations) were also constructed and tested for best fit among all other models. Finally, we also included a model predicting hummingbird species richness based on the log-transformed (base 10) total area in kilometers-squared of each forest fragment and, for illustrative purposes, compared it to a species-area relationship of hummingbirds from island ecosystems across the Americas using published literature [36][37][38][39] and data provided by the ebird online birding community ( [40]; Dataset S1). All models were compared using the corrected Akaike's information criterion (AICc) in R version 3.4.4. The calculation of AICc includes a correction factor that accounts for sample size and is preferable over AIC when sample size is small. We report top model coefficients and statistics for each figure in the supplemental tables.
To measure dimensions of hummingbird interaction specialization (i.e., niche partitioning and niche breadth) in different communities, we partially reconstructed hummingbird-plant interactions at 19 of our study locations from all captured hummingbirds by collecting pollen by dabbing their bills, lores, chin, and forehead with a piece of scotch tape in a standardized way during capture [14]. Pollen collected on the tape was placed on a glass slide for identification. We contracted an expert botanist to identify all pollen samples to morphospecies or the highest possible taxonomic resolution (henceforth called "pollen types") and used the number of times each pollen type was present on a hummingbird species (regardless of abundance within each pollen load) to assess the interactions of each individual. All pollen identifications were completed by comparing pollen shape, size, and diagnostic features against a photographed reference collection we created for this study, which was also used in previous studies [15]. We built one hummingbird-pollen type interaction network for each of the 19 sites and one for the entire agricultural ecosystem by pooling all sites into a single large network. To measure interaction specialization, we used a standardized index for calculating niche partitioning of a complex network called H 2 ' [41], also referred to as a measure of network specificity [42]. The H 2 ' metric is the standardized two-dimensional Shannon entropy for interactions rather than species [41,43]. A larger H 2 ' value indicates more niche partitioning, specifically via reciprocal specialization between hummingbirds and their plant partners [44]. Because almost all network metrics are correlated with species richness, sampling effort, and interaction distribution [45], all observed metric values were standardized against the distribution of expected values calculated from 1000 randomized null networks (the z-score). We calculated z-scores using the equation z-score = (x − x)/(σ x ) where x and σ x are the mean and standard deviation of the null distribution for metric x, respectively. The null networks were created using the Patefield algorithm [46] as in [47].
We also calculated hummingbird species specialization using the metric d', the species-level equivalent of the metric H 2 ' [41], as the metric d' is also a quantitative measure of species specialization, or specificity, [42]. Like H 2 ', a larger d' indicates greater species specialization. We calculated species specialization in two ways, first at the site level for each species (species_site d' z-score) to test for differences in niche breadth within the same species among different sites. Site-level species specialization was z-score standardized to compare among different networks and only included species that occurred at least once in all three habitat categories (14 species overall). Second, we calculated species specialization (species d') at the ecosystem level based on pooling all sites into a single large network. All network metrics and null models were calculated in the R package bipartite version 2.4 [45].
We gathered biophysical trait data to test for evidence that body size or bill morphology conferred a dependence on forest habitats or an advantage for colonizing coffee plantations. We took field measurement and measured hummingbird museum specimens at the Carnegie Museum of Natural History. Hummingbird species traits included standard measurements of bill curvature, bill length, body size, and wing length [15,26]. We also calculated the forest dependency of each species by creating a forest dependency index, which was the proportion of individuals captured in forest fragments and reserves combined after adjusting for sampling effort (Table S1).

Results
We banded 3205 individual hummingbirds at 19 sites and captured 24 of the 31 species of hummingbird ever recorded in the region [48]. Hummingbirds frequently captured included the Rufous-tailed Hummingbird (Amazilia tzacatl) and the Green Hermit (Phaethornis guy). We also captured 54 White-tipped Sicklebills (Eutoxeres aquila), a hummingbird with an extreme bill curve known to be a specialist [49]. There was no correlation between hummingbird species richness estimates and absolute forest fragment size in the agricultural ecosystem (mean ± SD = 11 ± 5 hummingbird species per site, regardless of forest fragment size or site identity as a coffee plantation: Figure S2), despite evidence of a species-area relationship based on hummingbird species richness on 76 islands across the Americas (Dataset S1; Table S2; Figure S2).
Hummingbird species richness estimates at 33 sites in the agricultural ecosystem changed in proportion to the amount of tree cover within 30 m of mist nets at each site. Moreover, the best model included a binary factor distinguishing coffee plantation versus forest sites, resulting in different y-intercepts in log space (R 2 adj = 0.44, p < 0.001, n = 33; Figure 1A, Table S3). Hummingbird abundance at 33 sites was also significantly related to the amount of tree cover within 30 m of the mist nets, again, with an important binary factor for sampling in a coffee plantation versus forest sites resulting in different y-intercepts in log space (R 2 adj = 0.16, p = 0.027, n = 33; Figure 1B; Table S4). Importantly, hummingbirds were twice as abundant in coffee plantations (mean abundance ± SD = 247 ± 61 individuals) compared to sites located in forest fragments and reserves (mean abundance ± SD = 138 ± 78 individuals). High abundances in coffee plantations were largely based on A. tzacatl abundances, Diversity 2020, 12, 126 5 of 14 a species with strong evidence that it colonizes coffee plantations in this region and elevation band. When we excluded data from A. tzacatl, hummingbird densities in coffee plantations (mean abundance ± SD = 90 ± 28 individuals) fell within the range of densities observed in forest fragments and reserves (mean abundance ± SD = 102 ± 74 individuals).
Diversity 2020, 12, x FOR PEER REVIEW 5 of 13 SD = 138 ± 78 individuals). High abundances in coffee plantations were largely based on A. tzacatl abundances, a species with strong evidence that it colonizes coffee plantations in this region and elevation band. When we excluded data from A. tzacatl, hummingbird densities in coffee plantations (mean abundance ± SD = 90 ± 28 individuals) fell within the range of densities observed in forest fragments and reserves (mean abundance ± SD = 102 ± 74 individuals). Hummingbird biodiversity changed in composition between forest sites and coffee plantations, suggesting the formation of a novel biological community. Evidence of a novel hummingbird community was detected by comparing Bray-Curtis dissimilarity coefficients between binary categories of coffee plantations versus forest fragments and reserves using permutational MANOVA tests (perMANOVA R 2 = 0.40, p < 0.01, F = 5.30). To visually demonstrate the novel hummingbird community, we plotted Bray-Crutis dissimilarity in a nonmetric multidimensional scaling plot ( Figure 2A). Because the composition of hummingbirds was driven by changes in species and abundances, we also plotted hummingbird effort-adjusted proportions of unique individuals captured for each species in order of their dependency on forest sites ( Figure 2B). Hummingbird biodiversity changed in composition between forest sites and coffee plantations, suggesting the formation of a novel biological community. Evidence of a novel hummingbird community was detected by comparing Bray-Curtis dissimilarity coefficients between binary categories of coffee plantations versus forest fragments and reserves using permutational MANOVA tests (perMANOVA R 2 = 0.40, p < 0.01, F = 5.30). To visually demonstrate the novel hummingbird community, we plotted Bray-Crutis dissimilarity in a nonmetric multidimensional scaling plot (Figure 2A). Because the composition of hummingbirds was driven by changes in species and abundances, we also plotted hummingbird effort-adjusted proportions of unique individuals captured for each species in order of their dependency on forest sites ( Figure 2B). categories of coffee plantations versus forest fragments and reserves using permutational MANOVA tests (perMANOVA R 2 = 0.40, p < 0.01, F = 5.30). To visually demonstrate the novel hummingbird community, we plotted Bray-Crutis dissimilarity in a nonmetric multidimensional scaling plot (Figure 2A). Because the composition of hummingbirds was driven by changes in species and abundances, we also plotted hummingbird effort-adjusted proportions of unique individuals captured for each species in order of their dependency on forest sites ( Figure 2B).

Figure 2.
Analysis of community composition suggested that hummingbirds colonizing coffee plantations formed a novel biological community in the agricultural ecosystem. A multidimensional scaling plot of hummingbird abundance-based biodiversity coefficients (Bray-Curtis Dissimilarity) illustrate differences and similarities among sites, where sites located closer together on the plot were more similar in hummingbird biodiversity and farther distances between sites indicate dissimilarity. Note how forest sites are more spread out in space, indicating that they host more diversity across sites from a larger pool of species (A). A plot of effort-adjusted proportions of unique captures of each hummingbird species suggests substantial overlap between communities in coffee plantations and in forest sites, with 17 of 24 species observed in both habitat types. Six species were so forest dependent that they were never observed in coffee plantations, and 4 of those species were rare (<10 total unique captures). Hummingbirds are labeled by species codes (see Table S1) and numbers below species codes indicate the total number of unique captures from 19 sites. Forest dependency rank is based on the effort-adjusted proportion of unique captures in forest fragments and reserves (B).
We collected 2866 pollen loads from 2416 unique hummingbirds and 20 species. In total, we were able to establish 3614 unique hummingbird-pollen type interactions and to create networks for 19 sites. We identified 170 pollen types (26 species, 27 at the genus level, 16 at the family level, and 101 morphospecies). On average there were 1.84 ± 1.06 pollen types in each pollen load. In Table S1 we report the percentage of sampled pollen loads with zero pollen for each species, like the Purple-crowned fairy (Heliothryx barroti), known for piercing flowers. There was no pollen on 82% of the sampled loads from that species. Hummingbirds are not the primary pollinators of coffee flowers and coffee pollen was detected in only 15 pollen loads.
We constructed hummingbird-pollen type networks for 19 sites where pollen loads were sampled (Figure 3 and Figures S3-S6) and for the entire ecosystem by pooling all sites into a single large network. All 6 coffee plantation sites and 5 of the 9 forest fragment sites were dominated by interactions between A. tzacatl and Erythrina spp. (Figure 3D-E, Figures S3-S6). The pollen type Erythrina spp. included the species Erythrina peoppigiana, a tree cultivated with coffee. Networks from sites in forest reserves and some forest fragments were more diverse than coffee plantations ( Figure 3A-C) and all 4 sites in forest reserves and 4 of 9 sites in forest fragments were dominated by interactions between P. guy and Heliconia tortuosa ( Figure 3F-G). Erythrina spp. pollen was found on 15 hummingbird species and was observed at all 19 sites. Pollen from Heliconia tortuosa was found on 12 hummingbird species and was observed at all 19 sites. The most specialized hummingbird in the agricultural ecosystem, E. aquila, was mostly captured in forest sites ( Figure 3H) and pollen load sampling confirmed a strong interaction with Centropogon granulosus, a flower with the same extreme curve as the bill of E. aquila ( Figure 3I). interactions between P. guy and Heliconia tortuosa ( Figure 3F-G). Erythrina spp. pollen was found on 15 hummingbird species and was observed at all 19 sites. Pollen from Heliconia tortuosa was found on 12 hummingbird species and was observed at all 19 sites. The most specialized hummingbird in the agricultural ecosystem, E. aquila, was mostly captured in forest sites ( Figure 3H) and pollen load sampling confirmed a strong interaction with Centropogon granulosus, a flower with the same extreme curve as the bill of E. aquila ( Figure 3I). Boxes on the left are hummingbird species; boxes on the right are pollen types. Height of the boxes represents relative abundance of hummingbirds and frequencies of pollen type in each network. Width of the connecting lines represents the number of times an interaction between a hummingbird species and a pollen type was observed. Symbols for each hummingbird and pollen type indicate their locations in the representative networks. Networks from sites located in coffee plantations were dominated by A. tzacatl (D). A. tzacatl interacted with 150 pollen types and mostly with Erythrina spp., whose frequency at each site was correlated with the abundance of A. tzacatl, meaning that as catches of A. tzacatl increased, so too did the number of observations of Erythrina spp., indicating frequent interactions between the two. (E). Networks from sites located in forest sites were dominated by P. guy, a hummingbird whose abundance increases proportionally with increasing tree cover within 30 m of mist nets at each site (F). P. guy interacted with 64 pollen types and mostly with Heliconia tortuosa, whose frequency at each site was correlated with the abundance of P. guy (G). Hummingbirds that tended to be more forest dependent also tended to be more specialized, especially E. aquila, whose rare observations were mostly in forest sites (H). E. aquila interacted with 11 pollen types and mostly with C. granulosus, whose frequency at each site was correlated with the abundance of E. aquila (I). All trend lines in subfigures (D-G and I) represent best fit models with a p-value ≤ 0.001.
Hummingbird niche partitioning, calculated using the network specialization index (H 2 ' z-score), also changed with the amount of tree cover within 30 m of the mist nets at each site (R 2 adj = 0.75, p < 0.001, n = 19; Figure 4A; Table S5). Site-level species specialization (species_site d' z-score) was calculated for 14 hummingbird species that occurred in all three habitat types, resulting in 169 site-level values. A generalized linear mixed effect model (GLMM) of site-level species specialization and the amount of tree cover within 30 m of the mist nets at each site was significant, suggesting that some hummingbirds are less specific in their interactions and relax their specialized relationships when they live and/or forage in coffee plantations ( Figure 4B; Table S6). Because the model predicted relatively small increases in species specialization with increasing tree cover, we assessed whether changes were a result of a type I error by comparing our data to null models. To create a null model, we randomly shuffled the data within species 5000 times and ran a GLMM for each random shuffle. From these randomly shuffled datasets, we found a significant relationship between tree cover and species specialization in less than 5% of the models, supporting our result that hummingbird specificity varies with tree cover within 30 m of mist nets at each site. relatively small increases in species specialization with increasing tree cover, we assessed whether changes were a result of a type I error by comparing our data to null models. To create a null model, we randomly shuffled the data within species 5000 times and ran a GLMM for each random shuffle. From these randomly shuffled datasets, we found a significant relationship between tree cover and species specialization in less than 5% of the models, supporting our result that hummingbird specificity varies with tree cover within 30 m of mist nets at each site.  Analysis of species-specific specialization in the ecosystem-wide network (species d', Table S1) suggested that more forest dependent hummingbirds tend to be more specialized, while less specialized species tend to be found in coffee plantations (R 2 adj = 0.23, p = 0.019, n = 20; Figure 4C; Table S7). In our relatively small analysis of 20 hummingbird species we found no support that traits of bill curvature, bill length, or body size conferred any advantage for using coffee plantations or predicting species specialization ( Figure S7A-C). Each species specialization value (species d') seemed to match with reported life history observations for each [49,50].
For example, E. aquila with its extreme bill curve was the most specialized hummingbird species (species d' = 0.82) with 2% of unique captures occurring in coffee plantations. The abundant P. guy was the third most specialized hummingbird species (species d' = 0.44) with 12% of unique captures occurring in coffee plantations. The most common bird in the agricultural ecosystem, A. tzacatl, was ranked 13 out of 20 in terms of specialization (species d' = 0.26) with 75% of unique captures occurring in coffee plantations. Other hummingbird species that preferred coffee plantations had similarly low species-specific specialization values (species d'). For example, both the Scaly-breasted Hummingbird (Phaeochroa cuvierii) and the Long-billed Starthroat (Heliomaster longirostris) represented about 90% of unique captures occurring in coffee plantations, and these species scored low specialization values (species d' = 0.13 and 0.23, respectively).

Discussion
Forest reserves and other high-quality forest habitats supported very complex and diverse hummingbird communities ( Figure 4A). Despite the reduced perception of human-dominated habitats to support biodiversity, lower-quality forest fragments and coffee plantations also supported multiple dimensions of hummingbird biodiversity, suggesting that there are conservation opportunities in the working landscapes outside of the irreplaceable reserves and protected lands. Our multiscale analysis suggests hummingbird species richness responds to tree cover at fine spatial scales (about 30 m), as it positively and logarithmically increased with tree cover both in forested habitats and coffee plantations. Second, a community analysis indicated that the hummingbird species found in coffee plantations overlap substantially with those in forest fragments, but hummingbird communities differ mainly in changes in species' abundances among habitats ( Figure 3B).
Our findings support previous studies that some hummingbirds increase in abundance after deforestation [18], such as A. tzacatl, which seems capable of colonizing coffee plantations. Despite some species "winning" in the wake of deforestation, hummingbird biodiversity followed patterns observed in other taxa in tropical agricultural landscapes [51], such as less diversity between sites in coffee plantations compared to forest habitats ( Figure 3B) and more evolutionarily unique species, such as E. aquila ( Figure 3H, Figure 4C), in forest habitats. We also report that 18 of the 23 hummingbird species observed in the agricultural landscape used coffee plantations to a greater or lesser degree ( Figure 3B). These results invite questions about the role of forest fragments and forested elements in coffee plantations in supporting or altering pollination dynamics in agricultural landscapes [15]. In addition, since hummingbirds are not the primary pollinators of coffee flowers (that role is filled by insect pollinators), it is not the coffee itself that is contributing to the creation of these novel hummingbird communities in agricultural habitat. Rather, it is likely the overall landscape context, including reduced canopy cover, different ambient temperatures and humidity, a different suite of floral resources, among other factors that may all determine which hummingbird species are successful in coffee habitat and which are not.
The hummingbird-plant interaction networks suggest that the novel hummingbird community we observed in coffee plantations was overall less specialized than the hummingbird community found in forest habitats, especially in forest reserves ( Figure 4). Hummingbird-plant networks in coffee plantations displayed greater interaction overlap, demonstrating that hummingbirds are not partitioning floral resources or forming reciprocally specialized interactions like they are inside forests ( Figure 4A). Rather, the hummingbirds found in coffee plantations were more opportunistic or random in their foraging patterns, and within-species specialization declined with decreasing forest cover. We found little support for the idea that physical hummingbird traits (such as bill length or curvature) predict species-level specialization or dependence on forested habitat ( Figure S6). The exception appears to be E. aquila, a species with a strongly curved bill, highly specialized interactions, and a strong dependence on forest habitat. Overall, the lack of relationship between morphological traits and interaction specialization may be problematic for studies that aim to infer hummingbird specialization based on morphology alone [15,19,26].
While we do not have direct measurements of floral resource availability at our sites, pollen type richness was similar between forests and coffee plantations, providing evidence that floral resources are available to hummingbirds found in coffee plantations. In particular, pollen from Erythrina species was a resource for 14 of the hummingbird species we captured, and Erythrina poeppigiana is a tree species maintained by farmers on coffee plantations to adjust microclimates and improve soil nutrients [52]. This is an example of an on-farm conservation initiative that benefits humans and hummingbirds alike. Evaluating the density of floral resources in coffee plantations may provide an opportunity to explore how resources and energy availability determine species colonization of new habitats, instead of quantifying habitat only in terms of total area [53,54].
Declines in within-species specialization (d ) along the deforestation gradient ( Figure 4B) indicate that hummingbirds can change their interaction behavior and become less specialized in different environmental contexts. Changes in foraging behavior of pollinators in response to other factors such as resource availability and elevation have also been recorded [11,20,23,27,47,[55][56][57][58], providing further evidence that hummingbird species are generally labile in their interaction partner choice under different circumstances. The differences we report in interaction niche partitioning and abundances in coffee plantations lead to questions about species evolution in human-made habitats [59] and continuation of pollination services to native plants [10,15,29]. These results highlight the need to study species interactions in both natural and human-made habitats in order to capture the full spectrum of interaction flexibility and vulnerability.
In many other regions of the tropics and the world, land-use change in the form of agricultural development is much more extreme, with industrialized, high-input monocultures dominating the landscape. Despite not including the extreme end of the agricultural management spectrum, our study still found significant changes in diversity, community composition, and interaction specialization in response to deforestation. It is likely that these effects are amplified in more intensive agricultural regions where remnant forests and clusters of trees have been removed, with even more potentially severe consequences for biodiversity maintenance and the sustained delivery of pollination services to native plants. Promisingly, flowering plants and trees planted by farmers for nutrient enrichment of the soil, noncommercial food resources, and microclimate amelioration are likely supporting hummingbird abundance in agricultural sites in our study region. In this vein, plant-pollinator networks can be used to identify habitats and resources used by vulnerable pollinators to inform targeted on-farm plantings and small-scale land sparing initiatives intended to benefit wildlife where reserves are infeasible and in tropical regions where biodiversity is high. While interaction networks show promise for informing conservation efforts, more empirical studies on species interactions across the spectrum of agricultural management types are needed to ensure that conservation strategies fully comprehend the vulnerability and flexibility of species interactions in an increasingly cultivated world.
Supplementary Materials: The following are available online at http://www.mdpi.com/1424-2818/12/4/126/s1, Figure S1: We sampled hummingbirds and their pollen loads over 5 years at 19 sites using mist nets. One site is not shown, a forest reserve site inside La Amistad Biosphere Reserve (4010 km 2 ) located~25 km Northeast of this map. We also show site locations where hummingbirds were sampled inside 14 forest fragment, reported by [19] who sampled hummingbirds using 10 mist nets operating under the same protocols and similar time periods as our sampling. Pollen loads were not collected by [19] and their data were only used for creating species-area relationships for hummingbirds ( Figure 1A) and modeling total hummingbird abundance ( Figure 1B). Figure S2: The number of hummingbird species relative to absolute forest fragment size or island size. Observed hummingbird species on 76 islands ranged from 0 to 18 species, with most islands hosting 0 to 4 species (Dataset S1). Dashed line represents a hummingbird species-area relationship for island ecosystem based on observed hummingbird species + 1 for true islands. There was no significant correlation between Chao hummingbird species estimates and absolute forest fragment size in the agricultural ecosystem. In the agricultural ecosystem, Chao hummingbird species estimates ranged from 4 to 20 species with a mean ± SD of 11 ± 5 hummingbird species per site (indicated by the solid line). Figure S3: Hummingbird-pollen type networks of 6 coffee plantation sites with proportion of tree cover within 30 m of mist nets at each site between 6% and 31%. Boxes on the left are hummingbird species; boxes on the right are pollen types. Height of the boxes represents relative abundance of hummingbirds and frequencies of pollen type in each network. Width of the connecting lines represents the number of times an interaction between a hummingbird species and a pollen type was observed. Relevant information and statistics are listed below networks, including raw and standardized network specialization values (H 2 ' and H 2 ' z-score). Figure S4: Hummingbird-pollen type networks of 4 forest fragment sites with proportion of tree cover within 30 m of mist nets at each site between 55% and 89%. Boxes on the left are hummingbird species; boxes on the right are pollen types. Height of the boxes represents relative abundance of hummingbirds and frequencies of pollen type in each network. Width of the connecting lines represents the number of times an interaction between a hummingbird species and a pollen type was observed. Relevant information and statistics are listed below networks, including raw and standardized network specialization values (H 2 ' and H 2 ' z-score). Figure S5: Hummingbird-pollen type networks of 5 forest fragment sites with proportion of tree cover within 30 m of mist nets at each site between 91% and 97%. Boxes on the left are hummingbird species; boxes on the right are pollen types. Height of the boxes represents relative abundance of hummingbirds and frequencies of pollen type in each network. Width of the connecting lines represents the number of times an interaction between a hummingbird species and a pollen type was observed. Relevant information and statistics are listed below networks, including raw and standardized network specialization values (H 2 ' and H 2 ' z-score). Figure S6: Hummingbird-pollen type networks of 4 forest reserve sites with proportion of tree cover within 30 m of mist nets at each site between 91% and 97%. Boxes on the left are hummingbird species; boxes on the right are pollen types. Height of the boxes represents relative abundance of hummingbirds and frequencies of pollen type in each network. Width of the connecting lines represents the number of times an interaction between a hummingbird species and a pollen type was observed. Relevant information and statistics are listed below networks, including raw and standardized network specialization values (H 2 ' and H 2 ' z-score). The far right hummingbird-pollen type network represents the forest reserve site located inside the huge La Amistad Biosphere Reserve (4010 km 2 ). Figure S7: Species-level specialization values based on pooling all sites into a single large network (species d') were not significantly correlated with bill curvature (A), bill length (B), or body mass (C). A. tzacatl (RTAH), P. guy (GREH), and E. aquila (WTSI) are labeled with symbols and images from Figure 4. Open circles in the plot represent 17 other hummingbird species with enough data to calculate species specialization. The correlation between bill curvature and species-level specialization produced significant p-values during analysis, but the removal of E. aquila from the model yielded an insignificant correlation. Table S1: Hummingbird species abundances and trait information, including species specialization index values (species d'), bill morphology, body size, wingbeats, forest dependency index values (percentage of individuals captured in forest fragments and reserves combined after adjusting for effort), and pollen load information. NA represents unavailable data. Table S2: Island ecosystem model estimates and structure for hummingbirds. This equation is a traditional species-area relationship that models the relationship between hummingbird species (S) and island area measured in kilometers-squared (A). Model was built using 76 islands from the Caribbean, California Channel Islands, Lake Gatún in Panama, and Gurí Lake in Venezuela. Log base 10 was used to match traditional species-area relationships. During analysis, hummingbird species richness values were modified by adding 1 to facilitate log transformations (+1) and the model includes a factor (−1) in its structure to back-transform species richness estimates. Table S3: Estimates and structure of top-ranked linear model representing the relationship between hummingbird species (S) and the proportion of tree cover within 30 m of mist nets at each site in the agricultural ecosystem (tree cover within 30 m). The proportion of tree cover measured at other spatial scales was also significantly related to bird species richness, but competitive models were discarded after comparing the corrected Akaike's information criterion (AICc). Natural log was used to demonstrate that species-area relationships are logarithmic in form. Table S4: Estimates and structure of top-ranked linear model representing the relationship between hummingbird abundance (A) and the proportion of tree cover within 30 m of mist nets at each site in the agricultural ecosystem (tree cover within 30 m). The proportion of tree cover measured at other spatial scales was also significantly related to bird species richness, but competitive models were discarded after comparing AICc. Natural log was used to demonstrate that abundance-area relationships are logarithmic in form. Table S5: Top-ranked linear model from of the relationship between tree cover and hummingbird-pollen type network standardized specialization index, a measure of niche partitioning (H 2 ' z-score). Table S6: Top-ranked generalized linear mixed-effects model of the relationship between tree cover and site-level species specialization index, a standardized measure of diet breadth of each species at each site (species_site d' z-score). Table S7: Estimates and statistics of the linear model between hummingbird specialization index (species d') and the forest dependency index (the effort-adjusted proportion of unique captures in forest fragments and reserves). Table S8: The coordinates (in UTM) and elevation (in meters) of each of the 19 study sites sampled for hummingbirds and hummingbird-plant interactions. The station code is a four letter code name of the sample site. Dataset S1: A comma separated value table listing the source, island identification in source, and reported hummingbird species richness for 61 islands from published literature [36][37][38][39] and 15 large (>500 km 2 ) Caribbean islands provided by the ebird online birding community, accessed on January 2019 [40].