Complementary Effects of Species Abundances and Ecological Neighborhood on the Occurrence of Fruit-Frugivore Interactions

Introduction

The structure and dynamics of ecological systems depend on plant-animal mutualistic interactions (e.g., pollination, seed dispersal), which are crucial processes underpinning community stability and ecosystem functioning (Okuyama and Holland, 2008; Zhang et al., 2011; Vázquez et al., 2015). As such, networks of mutualistic interactions are considered to be the architecture of biodiversity by contributing to its maintenance and stability and affecting coevolutionary processes (Bascompte and Jordano, 2007). Based on these observations, there has been growing interest both in describing the spatio-temporal patterns of pairwise interactions and in understanding the drivers of interaction occurrence that ultimately generate this architecture (Kaiser-Bunbury et al., 2014; Vizentin-Bugoni et al., 2014; González-Castro et al., 2015; Bartomeus et al., 2016).

Mutualistic interactions are fundamentally dynamic, such that there is a turnover of species and individuals in space and time (Carnicer et al., 2009; Laliberté and Tylianakis, 2010; Hagen et al., 2012). Despite this turnover, there may be some predictable patterns in the spatio-temporal distribution of pairwise interactions that are reflected in the emergence of areas/periods of time that are characterized by a higher abundance of interactions (Blendinger et al., 2015; Gilarranz et al., 2015). Identifying these spatial patterns could be valuable, as they may establish a template for the distribution of the ecological outcomes of pairwise interactions (e.g., higher plant recruitment derived from the activity of seed dispersers; Blendinger et al., 2008; Hampe et al., 2008). Indeed, the spatio-temporal patterns of pairwise interactions may condition relevant processes such as density-dependent mortality, genetic exchange, and ultimately the population dynamics of interacting species (Rodríguez-Rodríguez et al., 2015).

It is generally agreed that both neutral and niche-driven processes contribute to the occurrence of pairwise mutualistic interactions and network structure (Bartomeus et al., 2016; García, 2016; González-Varo and Traveset, 2016; Sazatornil et al., 2016). That is, whether or not two species interact can be thought of as the result of two sequential processes. First, there is a stochastic (neutral) process driven by the local abundances of interacting species, such that abundant species are more likely to encounter each other than are rare ones (Vázquez et al., 2007; Krishna et al., 2008). Secondly, there is a deterministic process driven by a matching of traits (e.g., phenology, morphology, behavior, etc.) that finally permits the interaction between co-occurring species (Maglianesi et al., 2014; Gonzalez and Loiselle, 2016; Donoso et al., 2017). Although both mechanisms are not mutually exclusive, neutral processes are usually expected from systems dominated by generalist species (e.g., González-Castro et al., 2015), whereas niche-driven processes are expected to be more relevant in those systems where high biodiversity involves a large trait variation among species (e.g., in tropical contexts, Bender et al., 2017). Unfortunately, a focus restricted to the dichotomy between these mechanisms ignores other drivers related to the occurrence or abundance of species other than the interacting pair in the community (e.g., context dependency driven by competition or facilitation by other species) (Perea et al., 2013; Martínez et al., 2014; Albrecht et al., 2015). These context dependencies, also termed neighborhood effects, may be inferred from the occurrence of third species comprising the spatial or temporal biotic context surrounding pairwise interactions. For example, in the case of interactions between fruiting plants and frugivorous animals, the spatial patterns of species occurrences have suggested that frugivore species may track each other (Saracco et al., 2005) or be attracted by large multispecies fruiting patches (Martínez et al., 2014). However, to our knowledge, no other empirical study has considered the presence of other co-occurring species to predict pairwise interaction frequencies (but see Kaiser-Bunbury et al., 2014; Poisot et al., 2015, for theoretical frameworks for doing so).

Here, we evaluate the patterns and drivers of fruit-frugivore interactions at a landscape scale over three consecutive years. We have two specific aims. First, to describe the spatial distribution of pairwise interactions and quantify their spatial consistency across years, together with a characterization of the assemblages based on the relative frequency of those interactions. This is a preliminary, but necessary step to the mechanistic understanding of pairwise interactions. Second, to estimate the simultaneous importance of neutral (i.e., abundances of interacting species) vs. deterministic (i.e., species identities) determinants, and other direct and indirect effects (i.e., neighborhood) on the occurrence and frequency of pairwise interactions across years using structural equation modeling. As a study system, we focus on an assemblage of three fleshy-fruited tree species and six frugivorous thrushes in a temperate forest of northern Spain. Previous studies in the same study system have demonstrated strong inter-annual variability in fruit production, bird abundances, and interaction diversity (García et al., 2013; Donoso et al., 2016), as well as the occurrence of indirect interactions between neighboring fruiting trees (Martínez et al., 2014). Amidst this variability, we aim to determine whether variability in the mechanisms driving pairwise interactions generates spatio-temporal patterns in the distribution of fruit-frugivore interactions. Although this study does not seek to evaluate how interaction patterns may influence demographic processes in populations, it will provide new insight into whether or not there is a spatial distribution of pairwise interactions, which mechanisms could generate these interactions, and how they might change over time.

Materials and Methods

Study System

This study focused on the mutualistic assemblage of fleshy-fruited trees and frugivorous birds in the temperate secondary forest of the Cantabrian Range (northern Iberian Peninsula). Due to anthropic pressure, this habitat is characterized by a high degree of fragmentation as well as low forest cover; specifically, it is predominantly secondary forest intermingled with much less frequent variable-sized fragments of primary forest, within a dominant non-forested matrix (García et al., 2005a). The fleshy-fruited trees selected for the study were the dominant species holly (Ilex aquifolium), hawthorn (Crataegus monogyna) and yew (Taxus baccata), which accounted for ca. 70% of tree cover in the studied forest (García and Martínez, 2012). All of these species ripen in autumn, with an overlapping period from September to November, but their fruits stay on the tree until mid-winter (January–February). Their fruits are sugar-rich red berries (arillated seeds in the case of yew) with similar morphology, coloring and size. Their main frugivores are thrushes (Martínez et al., 2008): blackbird (Turdus merula), song thrush (T. philomelos) and mistle thrush (T. viscivorus), which are resident species, and fieldfare (T. pilaris), redwing (T. iliacus), and ring-ouzel (T. torquatus), which are over-wintering species in northern Spain. Despite sharing a strong taxonomic affinity, these thrush species may differ in their response to landscape structure and habitat features (Martínez et al., 2008; Morales et al., 2013).

Study Site

The study area was located in the Sierra de Peña Mayor 43°18′00″N, 5°30′29″W, 1000 m a.s.l., Asturias (Spain). This area consists of a mountain range where secondary forests are juxtaposed with mature forest and a historically deforested matrix composed of heathland, meadows and rocky outcrops. We set up a rectangular plot of 400 × 440 m (17.6 ha) covering a secondary forest with a gradient of forest cover, from dense covered sectors to remnant trees, and this plot was then subdivided into 440 sampling cells of 20 × 20 m (see Supplementary Figure 1). Previous studies in the same area support the appropriateness of this spatial extent to represent the landscape patterns of fruit production, bird activity and forest cover (García and Martínez, 2012; Martínez and García, 2015a).

Monitoring Pairwise Interactions

From October to February of the three sampling seasons corresponding to 2008–2009, 2009–2010, and 2010–2011 (2008, 2009, 2010 hereafter), and from five vantage positions situated along the central axis of the plot (see Supplementary Figure 1), we observed the foraging behavior of birds during 90, 79, and 63 h (for 2008, 2009, and 2010, respectively). One-hour observation slots were equally allocated to the different vantage positions through the whole sampling season. For every frugivory event (a bird picking fruits from a tree) recorded during observation slots, we recorded the identities of the thrush and plant species and the number of fleshy fruits consumed per individual bird. We defined a pairwise interaction as one fruit of a given tree species removed by a given bird species (e.g. C. monogyna—T. merula; I. aquifolium—T. philomelos; etc.). We obtained the total pairwise interaction frequency per sampling cell based on the cumulative number of fruits of each tree species consumed by each bird species by the end of each sampling season. As different cells received different observation times (some cells were observed from several vantage points), we estimated pairwise interaction frequency per cell as the cumulative number of fruits consumed per 10 h of observation.

Fruit and Bird Abundances

In October of each sampling year, we mapped the position of all individual trees within the plot. The crop sizes of all fruiting individuals of the three species under study were visually estimated by means of a semi-logarithmic Fruit Abundance Index defined by six intervals (FAI, 0: no fruits; 1: 1–10 fruits; 2: 11–100; 3: 101–1,000; 4: 1,001–10,000; 5> 10,000; Saracco et al., 2005). Crop size was then extrapolated following an allometric equation (crop size = 1.77 × e^1.92FAI; R² = 0.80; n = 136 trees; Herrera et al., 2011). Each year, we incorporated this information into a GIS platform to obtain total fruit abundance per species per cell. In our system, fruiting of all individuals of the different species is synchronous and fruits ripen within 1–2 months. Therefore, we considered that a single sampling of fruit abundance at the beginning of the season would provide a suitable estimate of the spatial template of fruit resources available to frugivores (García et al., 2013; Martínez and García, 2015b). This template represents a static image of the maximum fruit availability at the beginning of the season, assuming that relative differences in fruit abundance between areas remain constant over the season. Although our use of this early-season measure could limit the analysis regarding the mechanisms driving late in the season interactions, we expect it to have a negligible effect on our global inference since most of the interactions were recorded early or mid in the season (75% of fruit consumptions where recorded by mid-December in 2008 and by late November in 2009 and 2010).

Bird counts were performed from the five vantage points by means of 1-h observation slots allocated equally across points (cumulative observation time of 103, 105, and 156 h, for 2008, 2009 and 2010, respectively), from October to February of each sampling year. Owing to lower bird detectability in some highly forested cells of the plot, we performed additional 10-min bird counts in 12 forest positions, with each one corresponding to the center of a group of four cells (see Supplementary Figure 1). In this case, the cumulative observation time from each point was 160, 110, and 195 min (for 2008, 2009, and 2010, respectively). The cumulative number of individual birds (thrushes) observed through the whole sampling season was estimated for each cell, and bird species abundance per cell was calculated as the number of birds per 10 h of observation, given that different cells received different observation efforts. These censuses led to an estimation of the activity of each frugivore species throughout the season for the different cells of the study plot, providing a measure of bird functional abundance. In those cases in which consecutive sightings of a given species could correspond to the same individual in the same cell or entering a given cell, we considered those records separated by at least 5 min as independent (see also García et al., 2013; Martínez and García, 2015a for details on bird count methodology).

Data Analysis

Data Spatial Pooling

For each year, we summed the data of fruit and bird abundance and pairwise interactions into groups of four adjacent cells, resulting in a total number of 110 different 40 × 40 m blocks for the whole study plot. This enlargement of the sampling units for analysis enabled us to work with a spatial extent large enough to detect the maximum values of species and pairwise interaction richness (i.e., the 6 species of birds, the 3 species of plants and, thus, the 18 potential pairwise interactions between them) (García and Martínez, 2012), as well as to represent the ecological neighborhood on a proper scale (see below and Martínez et al., 2014). In each block, we considered the combination of any fruiting species and any bird species that were both present (i.e., with non-zero abundance) to be potential pairwise interactions. Thus, for each block we assigned a pairwise interaction frequency of zero to any of these potential interactions that were not observed. On the other hand, we discarded from the analysis those blocks without birds and/or fruits recorded. Out of the total 110 blocks, interactions could occur in 80, 78, and 82 blocks in 2008, 2009, and 2010, respectively.

Quantifying Context Variables

We selected four variables to describe the surrounding environment of each pairwise interaction at the block scale, (i.e., the ecological neighborhood): (i) forest cover, estimated as the percentage of tree cover per block, estimated from a Geographical Information System (GIS, ArcGIS v9.3) of the plot based on 1:5,000–scale orthophotomap, incorporating a grid and a digitalized layer of tree canopy cover, including both forest patches and isolated individual trees; (ii) fruit patchiness, a measure of the spatial aggregation of fruit abundance based on a standardized clustering index (v), which quantifies the degree to which the abundance of fruits within each block contributes to the overall degree of fruit patchiness in the whole plot (from Spatial Analysis by Distance Indices; estimated by software Sadie Shell v.2.0 based on the sampled fruit abundances within cells; Perry et al., 2002). This variable has the advantage of reflecting, in relative terms, the total abundance of fruits as well as the contribution of each block to a total landscape-scale fruit aggregation; (iii) abundance of heterospecific fruits, the abundance of fruits from species besides that taking part in each pairwise interaction (iv) abundance of heterospecific birds, the abundance of birds from species other than that taking part in each pairwise interaction. Although there might be a temporal decoupling between the occurrence of each pairwise interaction and the surveys of fruit and bird abundances (in terms of bird activity), these estimates provided an initial template that is maintained (in relative terms) over the entire sampling season. Hence, the effects of heterospecific abundances on pairwise interactions would be mediated by a probability of encounter (e.g., the higher the abundance of heterospecific birds detected in one block through the whole season, the higher the probability of encounter between the birds of a given species visiting that block and these heterospecific birds).

Statistical Analysis

Identifying spatial patterns and determining their consistency through time

We were first interested in describing and quantifying the yearly spatial patterns of pairwise interactions. For this, we represented the percentage of pairwise interaction frequencies accounted for by each block within each year across the whole plot by means of a raster layer created with the raster package 2.5-8 (Hijmans, 2016). In order to assess the temporal consistency in the spatial distribution of pairwise interactions, we checked the correlation of the pairwise interaction frequencies across years, by means of the Kendall's coefficient of concordance (W) after 1,000 permutations (kendall.global function implemented in vegan R package 2.4-1; Oksanen et al., 2016). This coefficient measures on a 0 (no agreement) to 1 (perfect agreement) scale the degree of concordance among rank-ordered variables (pairwise interaction frequencies per year) for n objects or sites (110 blocks) (Legendre, 2005). We further performed an a posteriori test to assess the contribution of each individual year to the overall concordance (W per year) using the kendall.post function within the same R package. Finally, we also assessed the patterns of the spatial distribution of the total pairwise interaction frequency each year, by calculating an aggregation index (Ia) which represented a measure of the degree of spatial aggregation of pairwise interaction frequency at the plot scale (values of Ia = 1 representing random; Ia > 1 aggregated and Ia < 1 regular distribution patterns; from Spatial Analysis by Distance Indices; estimated by software Sadie Shell v.2.0; Perry et al., 2002). Since all the information was gathered in a spatially-explicit way, any non-random pattern in how interactions are distributed could help to elucidate which mechanisms were driving pairwise interactions.

Quantifying the mechanisms determining fruit-frugivore interactions

We employed a structural-equation-modeling approach to disentangle the effects of different mechanisms in determining the occurrence and intensity of pairwise interactions. We sought to represent three types of mechanisms through three different sets of variables, estimated for each pairwise interaction measured in each block each year. A first set of variables included the abundances of the bird and fruiting-tree species participating in the pairwise interaction, and were considered as the variables representing neutrality. A second set of variables accounted for the identity of the pairwise interactions (and that of species involved), as a measure of the deterministic effects of species trait-matching and behavior. For this, we firstly defined all the 18 different combinations of pairwise interactions considering the three tree species and the six bird species, generating a matrix of presence/absence data of the 18 pairwise interactions in each block each year. In order to represent the relative role of each pairwise interaction in the global trends of variability in the identity of interactions across the whole plot, we performed correspondence analyses (CA; Legendre and Legendre, 2012; see Supplementary Table 1 and Supplementary Figure 2) on the yearly interactions per block matrix. CA provided, for each pairwise interaction, scores for two principal dimensions that were used as continuous surrogates of the identity of pairwise interactions occurring each year across the study plot. CA also provided values of the relative contribution (proportion to inertia) of each pairwise interaction to the global yearly variability. We performed the CA using the ca R package (Nenadic and Greenacre, 2007) and the facto_summarize function from factoextra R package 1.0.3 (Kassambara and Mundt, 2016). Finally, the third set of variables represented the ecological neighborhood by means of forest cover, fruit patchiness, heterospecific bird abundance and heterospecific fruit abundance.

All eight of these variables were considered as potential predictors of the pairwise interaction frequency of different species pairs in each block of the study plot (excluding those blocks without potential interactions) in structural equation models (SEM; Shipley, 2009). We checked that all these variables were uncorrelated (R² < 0.4 in all cases). The SEM approach represents a form of path analysis that deals with complex multivariate relationships among a set of interrelated variables, and thus allowed us to consider the sum of causal direct and indirect interactions among variables (Grace, 2006). Path diagrams (the variables and effect pathways in Figure 2) were defined based on realistic causal relationships, using previous knowledge of our ecological system (e.g., Martínez and García, 2017). Specifically, we expected pairwise interactions to be shaped by abundances of potentially interacting species, by indirect effects of some factors on species abundances (e.g., heterospecific fruits, heterospecific birds, fruit patchiness and cover on abundance of interacting birds), and also by direct effects of these factors and species identities on pairwise interactions. As we were interested in comparing the same biological hypotheses across years, we performed the same SEM per year, with similar sample sizes across years (N = 536, N = 580; and N = 509 in 2008, 2009, and 2010, respectively; where N is the total number of pairwise interactions considering both realized interactions based on fruit consumption and potential pairwise interactions corresponding to the combination of fruit and bird species whose abundances were recorded in the block but whose interaction was not detected). Our SEMs were represented by a list of models considering pairwise interaction frequency and interacting-bird-species abundances as the two response variables. Compared with traditional variance-covariance based SEM, the piecewise SEM procedure allowed us to fit generalized linear models, considering non-Gaussian error distributions in response variables (Lefcheck, 2016). Namely, both the abundance of interacting birds and the pairwise interaction frequency were over-dispersed count outcome variables with many zero values, so we built models considering negative binomial distributions (by means of the glm.nb function from the MASS R package). Prior to analysis, all the explanatory variables were standardized and the goodness-of-fit of each model was also confirmed with a chi-square test based on the residual deviance and degrees of freedom. We further extracted all the standardized path coefficients for the structural equation model. These coefficients represent the estimate of each causal path between each response variable (i.e., pairwise interaction frequency and interacting bird species abundances) and each of the predictor variables. They were calculated by using the sem.coefs function from piecewiseSEM package. We further computed indirect effects of each predictor on pairwise interaction frequency as the sum of products of the coefficients along all of the possible routes from each predictor to the response variable.

Lastly, due to the spatial configuration of our sampling framework (i.e., 110 adjacent blocks), the estimation of the effects of interacting species abundances, interaction identity and neighborhood context on pairwise interactions, as well as on the abundance of interacting bird species, might have been affected by potential spatial autocorrelation. Thus, for each sampling year we also estimated spatial autocorrelation of the residuals of the list of models representing our SEM by means of Moran's I, using the correlg function from the ncf R package 1.1-7 (Bjornstad, 2016).

Results

Spatio-Temporal Patterns in Fruit-Frugivore Interactions

Field sampling revealed that the potential for pairwise fruit-frugivore interactions (i.e., expected from fruiting tree and bird presence) was widespread across the whole plot all years, with three quarters of blocks of the sampling plot harboring some potential interactions in each year. Nevertheless, the spatial distribution of the pairwise interaction frequencies (no. fruits consumed/10 h) per block changed across years (Figure 1, upper panels). The SADIE-based aggregation index revealed an overall aggregated pattern of fruit-frugivore interactions for 2008 (Ia = 1.50; p = 0.006) and 2010 (Ia = 1.33; p = 0.04). Nonetheless, pairwise interactions were aggregated in different areas of the study plot each year: they were mainly biased toward forest areas in 2008 but toward open deforested areas in 2010. It is worth noticing that the maximum value of aggregation also changed across years. For example, although interactions in 2008 showed an aggregated pattern in forest areas, there was a single block with an extremely high aggregation value. In 2009, the aggregation index suggested a less aggregated, more regular distribution across the whole study plot (Ia = 0.86; p = 0.79). Kendall's coefficient of concordance (W) of pairwise interaction frequencies per block across years revealed significant, but moderate, consistency in the spatial distribution of interactions during the studied period (W = 0.45; p = 0.006). This concordance probably resulted from the consistency in those areas encompassing zero or just a few interactions during the three sampling years, rather than from consistency in the distribution of the scarcer interaction hotspots (i.e., concordance disappeared when removing blocks encompassing zero interactions during the three sampling years; W = 0.32; p = 0.54). Moreover, the a posteriori test for the contribution of individual years to the overall concordance showed that the spatial distribution of the interactions in 2010 was less concordant with the other 2 years (see Supplementary Table 2).

Figure 1. In the upper panels, we show the spatial distribution of the sum of pairwise interaction frequencies accounted for by each of the 110 block (fruits consumed/10 h of observation; expressed as a percentage scale) in the study plot (blocks depicted by black grid; original cells depicted by white grid) during three consecutive years: (A) 2008, (B) 2009, and (C) 2010. The color scales are shown below each panel, but note that the values change between years. Blocks without foreground color represent areas where no interaction could occur as there were no fruiting trees and/or birds recorded during that particular year. In the lower panels, we show the distribution of frequencies of pairwise interactions corresponding to different birds, in a decreasing order (bird species abbreviation codes in the y axis; Tuil, T. iliacus; Tume, T. merula; Tuph, T. philomelos; Tupi, T. pilaris; Tuto, T. torquatus; Tuvi, T. viscivorus) and tree species (green hues, I. aquifolium; red, C. monogyna; blue, T. baccata) for different years in the study plot.

Likewise, we found interannual variability in the structure of the assemblages based on the relative frequency of each pairwise interaction (Figure 1, lower panels). In 2008 and 2009, I. aquifolium was the most consumed plant species, mainly by T. iliacus and T. merula, which together accounted for more than half of the total fruit consumption. In 2010, the frequencies of the different fruit-frugivore interactions were more evenly distributed, although they were predominantly made up by C. monogyna, and T. baccata consumed by T. viscivorus and T. merula.

Mechanisms Determining Fruit-Frugivore Interactions

Our evaluation of the influence of species abundances, interaction identities and the ecological neighborhood revealed an overall consistency across years in the main predictors of pairwise interactions (Figure 2). While the structural equation models included no statistically significant effect of interaction identities, the positive effect of both interacting plant and bird species abundances on pairwise interaction frequencies was consistent across years. The abundance of interacting bird species was also positively related to two neighborhood variables: forest cover and heterospecific bird abundance. In 2008 and 2010, fruit patchiness negatively affected interacting bird species abundance, and heterospecific bird abundance was negatively related with pairwise interaction frequencies. Finally, pairwise interaction frequency was only affected by a direct negative effect of heterospecific fruit abundance in 2008. Overall, we found that pairwise interaction frequency depended on the interacting species' abundances and the negative influence of heterospecific bird abundances, as well as other variables representing neighborhood context (i.e., cover or fruit patchiness), whose effects were only indirectly mediated by the abundance of interacting birds. The identity of interactions had no effect on pairwise interaction frequencies.

Figure 2. Panels on the left illustrate path diagrams representing how ecological neighborhood influences the abundance of the interacting bird species and the pairwise interaction frequency (based on fruit consumption data) during 2008, 2009, and 2010. The effects of both interacting bird and fruit abundances and that of species identities (based on the two dimensions of a Correspondence Analysis on matrices of presence/absence of pairwise interactions across study blocks) on pairwise interaction frequency are also shown. Solid arrows represent statistically significant relationships (p < 0.05); dashed lines are non-significant. Black and gray arrows represent positive and negative relationships, respectively. Values of standardized path coefficients are shown next to pathways for the significant relationships. Note that the thickness of significant paths has been scaled based on the magnitude of the standardized regression coefficient and can differ between years. Panels on the right show bar plots with the sign and magnitude of standardized direct (filled bars) and indirect (striped bars) effects of each variable predicting pairwise interaction frequency. Dashed lines separate the three sets of variables representing the ecological neighborhood (on top), deterministic effects (in the middle), and neutral processes (on the bottom).

Moran's I correlograms showed that both raw data and residuals of the model predicting pairwise interaction frequency were barely autocorrelated in space. Conversely, raw data for the abundance of interacting bird species presented significant spatial autocorrelation at short distances (up to 80 m). Nevertheless, this autocorrelation decreased strongly when considering model residuals, suggesting that our model incorporated the main processes behind the spatial patterns of interacting bird species abundances (McIntire and Fajardo, 2009; see Supplementary Figure 3).

Discussion

Here we have explored patterns and the underlying processes of pairwise interactions in a small assemblage of fleshy-fruited trees and frugivorous thrushes in a landscape of the Cantabrian range over three consecutive years. Our results illustrate how the distribution of pairwise interactions across species of fruiting plants and frugivores, as well as the large-scale spatial distribution of these interactions, changed across years. However, despite these changes, we found temporal consistency in the mechanisms responsible for the occurrence of pairwise interactions. In this sense, we used an integrative approach to disentangle the role of neutral (i.e., abundance-based) vs. deterministic (i.e., species-identity-based) drivers of interactions. Our study overcomes the limitations of classical views of fruit resource tracking (e.g., Blendinger et al., 2015) as it explains frugivory patterns from fruit and frugivore abundances by taking explicitly into account which bird interacts with which plant. More importantly, we have extended upon previous studies that used the abundance and traits of the interacting species (González-Castro et al., 2015; Bartomeus et al., 2016) by considering the effects of ecological neighborhood (including the indirect effects mediated by the presence of other species co-occurring with those of the interacting pair). Overall we found that the assemblage of pairwise interactions depended on the neutral effects of interacting species abundances (the more abundant the species of fruit and frugivore, the more frequent their pairwise interaction) and on the local context of ecological neighborhoods.

Interannual Variability in the Spatial Distribution of Pairwise Interactions

In this study, we found no evidence of a strong temporal concordance in the spatial distribution of interactions across years. The occurrence of areas with a high frequency of interactions is expected from the common non-random distribution of fruit resources and the non-random habitat use of frugivores (Jordano, 2000). Moreover, as fruiting trees are sessile organisms and frugivores frequently show preferences among the available fruiting plants (Carlo et al., 2007), temporal consistency in the location of high-frequency interactions could also be expected. Nevertheless, this does not seem to be the case in our study system. In fact, the interannual variability in the distribution of pairwise interactions seemed to be more strongly related to fruiting plants than to frugivorous thrushes. Specifically, the three species of trees under study are non-overlapping masting species with different standing distributions across the landscape. Therefore, different fruiting landscapes were expected in each year, generated by a dominant masting species and a different distribution of fruiting patches (Herrera and García, 2009; García et al., 2013; Donoso et al., 2016). This was certainly apparent in our study system, where fruits of I. aquifolium were more abundant and biased toward forest areas in 2008 and 2009, whereas the fruit crop in 2010 was dominated by C. monogyna, which was mainly present in open deforested areas (see Supplementary Figure 4). This variability in fruit-crop composition provoked by masting events of the main fruiting plants I. aquifolium and C. monogyna was also reflected in the different composition of pairwise interaction communities (Figure 1).

On the other hand, birds tracking those different fruiting templates each year would also lead to a different spatial distribution of interactions. From the perspective of our whole bird community, frugivory activity seemed to depend on forest cover and fruit availability at the landscape scale (Martínez et al., 2014; Martínez and García, 2015a). Furthermore, at a species level, the different thrushes are also known to respond differentially to landscape features (García et al., 2013; Morales et al., 2013; Rodríguez-Pérez et al., 2014) so the spatial patterns of interactions might also be affected by the different species identities composing pairwise interactions each year (Figure 1; see also Blendinger, 2017, where spatial changes in the outcome of interactions involving different species of thrushes were shown). For instance, while T. iliacus, a species typically restricted to large forest patches, was much more frequent in species interactions in 2008, T. viscivorus, a species more prone to explore the whole landscape, accounted for a large number of interactions in 2010. Therefore, the fact that generalist species could reorganize their interactions by replacing their partners suggested that species in our generalized plant-frugivore assemblage were efficient when tracking fruit resources and showed strong context dependency despite the interannual spatial variability (García et al., 2013; Chamberlain et al., 2014). Since seed removal has been shown to not be strongly affected by possible shifts in the frugivory community (Farwig et al., 2017), temporal inconsistency in the spatial patterns of interaction might have long-term effects on seed-dispersal processes, especially in such degraded and fragmented systems. Although denser seed rain is expected in areas concentrating pairwise interactions, as suggested by the typically higher density of dispersed seeds under fruiting canopies (García et al., 2005b; Herrera and García, 2010; Carlo et al., 2013), this temporal and spatial dynamism might be relevant for diluting density-dependent seed mortality, favoring genetic exchange, interspecific seed transfer (García et al., 2007) and promoting seed deposition and forest regeneration in open areas (Martínez and García, 2015b). Nevertheless, further studies are necessary to explicitly quantify the long-term demographic consequences of spatio-temporal patterns in pairwise fruit-frugivore interactions.

Consistency in the Mechanisms Determining Pairwise Interactions

We found that the abundance of plant and bird species was always strongly positively associated with the occurrence of pairwise interactions. Our measures of interacting species abundances were independent of those of pairwise interactions, so no spurious analytical correlation is expected (González-Castro et al., 2015). This positive effect of encounter probability was not consistent with that presented by Olito and Fox (2015) for plant-pollinator interactions, but agrees with other studies showing neutrality as an important predictor of species assemblages, together with niche processes following the continuum theory (Gravel et al., 2006; Krishna et al., 2008; Vázquez et al., 2009; Kaiser-Bunbury et al., 2014). Most of these studies (mainly on plant-pollinator interactions) disentangled the mechanisms behind the macroscopic, global features of interaction networks (Trøjelsgaard and Olesen, 2016), but failed when seeking to predict pairwise interactions (Vázquez et al., 2009; Vizentin-Bugoni et al., 2014; Gilarranz et al., 2015; Sazatornil et al., 2016). Here we found it to work on the individual interaction scale.

In the context of fruit-frugivore interactions, González-Castro et al. (2015) suggested that seed-dispersal interactions may be predicted by the combined effects of trait matching and species abundances. Hence it can be expected that strong trait matching between species would lead to a better prediction of pairwise interaction frequencies by species identities. Nevertheless, we did not find such effect in our study, given that the community-scale interaction frequencies were unaffected by variation in species identities across the landscape. Although some previous studies have suggested some kind of species-specific preferences derived from phenological matching (e.g., T. iliacus feeding on I. aquifolium or T. philomelos on T. baccata; Martínez et al., 2008), our results suggested that no strong trait matching or behavioral preferences were relevant when estimated in a species assemblage context. A stronger role of neutral mechanisms in predicting pairwise interactions was expected from our system, a typical temperate environment transformed by post-glacial human activity, usually composed of few and mostly generalist species, and in which all the bird species under study are congeners (i.e., they might exhibit less trait variability than those presented by González-Castro et al., 2015). However, trait matching could be expected to be a pivotal mechanism in those cases in which trait variation among species is larger, and where coevolutionary or long-term processes of ecological fitting have been operating (e.g., in tropical contexts; Dehling et al., 2014; Bender et al., 2017; Donoso et al., 2017). Nontheless, incorporating interaction identity in our model was essential since it would have revealed any possible preference between species.

In this study, we have also extended upon previous research by considering the role of variables representing the ecological neighborhood in determining pairwise interactions. This consideration mainly revealed that other co-occurring species could also play a role in predicting interactions (specifically by provoking relationships of competition between species). These neighborhood effects could be expanded if considering indirect effects derived from the constant and positive effects of both forest cover and heterospecific bird abundance on interacting bird abundance (facilitation between bird species). This suggested that the individuals of a given bird species tended to concentrate in areas of high forest cover but also in areas visited by other species. Although our analysis did not distinguish whether the presence of the different bird species in a given sampling block was simultaneous or consecutive, we argue these indirect effects to be mostly related with interspecific facilitation in habitat use (Saracco et al., 2005; Sridhar et al., 2009), as suggested by the frequent detection of mixed-species flocks in field observations. Surprisingly, in two of the three study years, we also detected a significant, negative direct effect of heterospecific bird abundance on the frequency of pairwise interactions. This means that the interactions involving a given bird species were rarer where other bird species were abundant, presumably due to resource competition. Hence, the positive relationship between bird species that we interpreted when birds are using the habitat (i.e., tracking each other across the landscape for arriving at the same sites) seemed to become negative (i.e., competition for fruit resources) once several bird species had met in the same patch. Fruit patch defense and winter territorialism, involving mostly large-bodied species like T. viscivorus (Snow and Snow, 1984; Skórka and Wójcik, 2005), would explain these competitive effects (which, in fact, were stronger in 2010 when paired interactions were dominated by mistle thrush; Figure 1). Finally, we also detected a negative effect of fruit patchiness on the abundance of interacting bird species, indicating that fewer interacting birds concentrated in those blocks with many fruits of all species that contributed the most to fruit patchiness. As fruit patchiness represents the contribution to global spatial aggregation of fruit crops, high values of fruit patchiness correspond to points within large fruiting areas, where frugivorous birds are probably widespread over larger extents leading to lower frugivore densities. This can help reconcile the negative effect of fruit patchiness on interacting bird species abundance. Furthermore, these high fruit patchiness points are also present in protective areas, where birds usually relax their fruit-tracking behavior, leading to unexpected null or even negative relationships between fruit availability and bird abundance (Martínez and García, 2015a). Overall, these results provide empirical support to the idea (Poisot et al., 2015) that the occurrence of a pairwise interaction is not only contingent on the presence of the species involved and their abundance, but also upon the neighborhood context and the effect of other species in the community that may influence local bird distribution (Blendinger et al., 2012). Likewise, other indirect neighborhood effects (e.g., forest cover, fruit patchiness or heterospecific bird abundances) mediated by the abundance of interacting birds may also play a role in predicting pairwise interactions.

Concluding Remarks

This study responds to the need to examine the variability of mutualistic relationships at the fine scale of pairwise interactions, while identifying the main forces (in our case species abundances but also environmental context) promoting such variability (Trøjelsgaard and Olesen, 2016). Despite showing different spatial patterns each year, we found that pairwise interactions were regulated by consistent mechanisms. Even though all species in the fruit-frugivore assemblage co-occurred in the landscape during the three sampling years, changing species abundances mainly explained the differences in fruit-frugivore interactions among years. Thus, we expect the derived ecosystem function provided by frugivorous birds, i.e., seed dispersal, to be determined by the spatial patterns of species abundances. Moreover, the use of a fine-grained spatio-temporally sampled data, together with an integrative approach provided by structural equation modeling, contributed to a more detailed mechanistic understanding of other extrinsic factors promoting, directly or indirectly, variability in the microstructure of our plant-frugivore assemblage. Therefore, we encourage future studies to incorporate such interactions that could also play a role through relationships of facilitation or/and competition among species, in order to improve the ability to predict the occurrence of pairwise interactions.

Author Contributions

Author sequence reflects decreasing order of contribution. DG and ID designed the study. DM and DG designed and conducted fieldwork. ID performed the analysis with inputs from DG, JT, and DS. ID wrote the first draft of the manuscript and all authors contributed substantially to revisions.

Funding

This research was supported by the Spanish National Program of Research and Development (CGL2015-68963-C2-2-R, CGL2011-28430 and PRI-AIBNZ2011-0863 grants to DG) and the FPI Program—European Social Fund (BES2009-25093 grant to DM, BES2012-052863 and mobility EEBB-I-14-08279 grants to ID). JT and DS were funded by Rutherford Discovery Fellowships, and DS by a Marsden Fund Fast-Start Grant (UOC-1101), all of which were administered by the Royal Society of New Zealand.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank J. Rodríguez, C. Guardado and R. Pérez for technical support during fieldwork. We are grateful to Jonathan S. Lefcheck and J. Rodriguez-Pérez for their statistical advice and to Soumya Prasad and Tomás A. Carlo for the opportunity to contribute to this special issue. Field work was conducted with the permission of the Wildlife Service of Asturias.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2017.00133/full#supplementary-material

References