Zooming into plant-flower visitor networks: an individual trait-based approach

Study system

We conducted fieldwork in a 0.52 ha flower-rich meadow (UK National Vegetation Classification MG5–unimproved old hay meadow) located at the Roding Valley Nature Reserve (Essex, UK [51°38′04.1″N 0°04′12.6″E]) (Fig. S1). The surrounding landscape comprised similar meadows with associated hedgerows and trees. Due to the labour-intensive effort required to sample traits for the interacting partners of the entire flower visitor community, we focused on bees, which stand out among insects as the world’s primary pollinators in most ecosystems (Neff & Simpson, 1993; Winfree, 2010). Our study site supported ten species (four families) of flowering plants and 28 species (six families) of bee (Tables S1, S2).

Sampling procedure

We conducted surveys over 5 days from 29 June to 5 July 2011, and 4 days from 13 June to 28th June 2012, using two observers to record interactions between bees and flowers along the transects. These short periods represent the peak period for the meadow which is cut for hay in July, thus removing all flowers. The transects were positioned to maximise coverage of all possible interactions in the meadow. In 2011 this resulted in an approximately ‘L’-shaped transect, with an angle of 105° between the two sections (Fig. S1A) and a total length of 135 m. To increase coverage of the meadow in 2012, we changed the transect into four separate sections (75, 53, 31, and 17 m) perpendicular to the length of the field with a total length of 176 m (Fig. S1B). Interactions were observed up to 1 m either side of the transects, meaning an effective coverage of 270 m² (5.2% of the total meadow area) in 2011 and 352 m² (6.8% of the total meadow) in 2012. To allow for diurnal turnover in interactions, the transects were walked twice per day in opposite directions between 10:00–13:00 and 12:00–16:00 in 2011. Due to increased duration of transect walks, and collection of other data, transects in 2012 were only walked once per day in opposite directions on consecutive days. Surveys lasted approximately 2 h per transect walk in 2011, with one transect conducted on the first day and two transects on all other days (approximately 15 h observations). In 2012, surveys lasted between 2 and 4 h. Total survey time across the two years was approximately 28 h. Surveys in both years were conducted during fine weather (the weather did not change during the course of any survey and no transects needed to be abandoned).

On each transect survey, one observer noted the type of interaction (nectar and/or pollen collection, or nectar robbing) for bees visiting flowers. Once the interaction was observed, bees were captured and the second observer measured the maximum floral display size of the flower being visited. We defined floral display size as the width at the widest point of a single flower, or in the case of Asteraceae and Apiaceae, a single inflorescence. We placed bees that could be identified in the field (e.g., Bombus spp.) in a bee-marking cage, positioned so that their bodies were flat and their tegulae were visible, and took digital photographs of them against a scale to allow the intertegular distance to be measured using image analysis software (ImageJ v1.46r). Prior to their release at the same location as they were captured, we marked all bees with a dot of non-toxic marker paint (Posca PC-5M, Uni, Japan) on the dorsal surface of the thorax to avoid any re-recording of the same individual. For any bees requiring laboratory identification, we used individual collection vials charged with ethyl acetate, and took a digital photograph of the pinned specimen to measure intertegular distance and proboscis length using the same image analysis software as above.

Both proboscis length and nectar holder depth are considered as pivotal functional traits in plant–pollinator interactions (Harder, 1983; Stang, Klinkhamer & Van der Meijden, 2006a; Stang, Klinkhamer & Van der Meijden, 2006b; Stang et al., 2009). Therefore, we used collected specimens of each bee species encountered (n = 1 to 29 depending upon encounter rate per species) to define the relationship between intertegular distance and proboscis length to enable us to use intertegular distance measured in the field to predict the proboscis length of bees from all observed interactions. We also collected a subsample of open flowers from each flowering plant species (5–10 flowers per plant, 1–5 plants per species; n = 5 to 25) to define the relationship between maximum floral display size and nectar holder depth per plant species. Regressions of nectar holder depth against maximum floral display size for each species were used to predict the nectar holder depth of flowers whose floral display size was measured during the transects. For species with nectar holder depths that were effectively too short to measure, we assigned a value of 1 mm for the nectar holder depth of all individuals (four out of 10 species). Due to the high abundance of certain bee species, such as honeybees (Apis mellifera) and bumblebees (B. lapidarius), and the time-consuming nature of trait measurements, individual measurements were only taken for the first 50 interactions observed in these cases. This allowed a wider range of unique individual interactions to be characterised.

Construction of flower-visitation networks

We used the frequency of flower visits by bees as the interaction weight to construct two different types of quantitative bipartite networks: (i) a traditional species-based network, and (ii) our novel individual trait-based network. For the trait-based network, we used the functional size measurements obtained from individual plants and bees. Predicted nectar holder depth (NHD; Fig. S2A) was used as the functional size trait to characterize our focal flowering plant community. We selected intertegular distance (ITD) as our measure of functional size in bees as it was highly correlated with proboscis length (Fig. S2B) and it captures other variables such as flight range (Greenleaf et al., 2007). This high-resolution individual data on trait variation (NHD and ITD) was used to perform cluster analyses and objectively allocate all interacting individuals into more realistic trait groups, which were then treated as nodes in the subsequent networks.

All analyses were performed in R 3.2.4 (R Development Core Team, 2015). As both NHD and ITD are continuous variables, we computed an agglomerative hierarchical clustering with the function agnes available in the cluster package (Maechler et al., 2016), using Euclidean distances for calculating dissimilarities between individuals. Agglomerative clustering starts with each individual contributing to the cluster analysis being treated separately, and then joins individuals into clusters based on the distance metric. Firstly, in order to compare the traditional species-based network with our trait-based network more directly, we ran the analysis while setting the number of clusters in the trait-based network to match the number of plant and bee species (following Woodward et al., 2010b) and remove the effect of network size on network properties. We refer to this hereafter as ‘constrained’ functional size-based network. Secondly, we explored the optimal number of clusters in which individual bees and flowers of the community group according to their trait variation (and irrespective of species), and constructed the resulting ‘unconstrained’ functional size-based network. To compute these cluster analyses, we used the NbClust package (Charrad et al., 2014), which provides up to 30 indices for determining the optimal number of clusters in the dataset. This approach offers the best clustering scheme by varying all combinations of number of clusters, distance measures, and aggregation methods. We selected the “average” aggregation method, which uses the average pairwise distance between all pairs of individuals in the different clusters as the measure of distance (Sokal & Michener, 1958). Network figures were drawn using the bipartite package (Dormann, Gruber & Fründ, 2008).

Network parameters and data analysis

We used four commonly used quantitative network parameters (weighted connectance, weighted nestedness, interaction evenness and degree of complementary specialization) and two node-level parameters (specialization and strength) to describe and compare the structure of the networks constructed:

Weighted connectance (C_q) (e.g., Kaiser-Bunbury et al., 2011): gives individual weight to each node based on their total interaction frequency, in contrast to qualitative connectance (the fraction of all possible links that are realized in a network), and, therefore, better captures the functional importance of a species (or functional group) in the community.

Weighted nestedness (WNODF) (Almeida-Neto & Ulrich, 2011): uses quantitative data to give a measure of the degree of hierarchy in the organization of the interactions. In our case this relates to size structuring of functional size based network. In a nested network, nodes with fewer interactions (specialists) are mostly linked with a subset of nodes linked to the most connected ones (generalists). Nestedness ranges from zero (not nested) to 100 (highly nested).

Interaction evenness (IE) (Tylianakis, Tscharntke & Lewis, 2007): is based on Shannon’s evenness and measures the uniformity of interactions among nodes in a network. IE ranges from zero (complete unevenness in the distribution of interaction frequencies) to one (complete uniformity).

Degree of complementarity specialization (H₂′) (Blüthgen, Menzel & Blüthgen, 2006): is a measure of specialization that depicts how much the interactions of each node differ from each other in the network. H₂′ ranges between zero (no specialization) and one (complete specialization).

Specialization (d′) (Blüthgen, Menzel & Blüthgen, 2006) for plants (d′_p) and bees (d′_b): gives levels of specialization of each species (or functional trait node in our case). It accounts for the available resources provided by the interaction partners, so a pollinator species (or group of species included in a node) that visits resources proportionally to the total number of interactions of the species (functional size nodes) it interacts with is considered generalized, while a species that visits rare resources disproportionately is considered specialized. The index ranges from zero (highly opportunistic) to one (highly selective).

Strength of a bee species (or node) (st′_b) is a measure of the importance of a pollinator from the perspective of the flowering plant community (Bascompte, Jordano & Olesen, 2006). It is the sum of dependencies of the plants relying on that particular pollinator. In the same way, the strength of a plant species (or species included in a node) (st′_p) is the sum of dependencies of the pollinators relying on that given plant species. In the case of the bees, for instance, it is calculated as the relative frequency of a bee species (or bee size cluster) on a particular plant species (or plant size cluster), i.e., the number of interactions between pollinator nodes j and plant nodes i divided by number of visits of all pollinator nodes to plant node i.

Network parameters were calculated using the bipartite package (Dormann et al., 2009) run in R 3.2.4. (R Development Core Team, 2015), testing their significance against 1,000 networks generated by the null model r2dtable (function ‘nullmodel’ in bipartite) based on the Patefield algorithm (Patefield, 1981), and using a z-score test. We also tested the significance for nestedness in our networks using three different null models (CRT, Conserve Row Totals; CCT, Conserve Column Totals; and RCTA, Row Column Total Average), implemented in FALCON (Beckett, Boulton & Williams, 2014). We used the adaptive ensemble method to reduce undersampling of the null distribution and to optimize the minimal ensemble size sufficiently to give robust statistics.

Finally, we also computed centrality to evaluate the node’s relative importance to the structure of the network and identify key species and functional size ranges in the community (Martín González, Dalsgaard & Olesen, 2010; Gómez & Perfectti, 2012; Mello et al., 2015). Centrality indicates how well connected a node is to the rest of the nodes in the network; central nodes can trigger a rapid breakdown of the network structure if they are selectively removed (e.g., Memmott, Waser & Price, 2004). We estimated centrality using three metrics (see Mello et al., 2015 for details): (i) normalized degree centrality, i.e., the proportion of different partners a certain node interacts with in relation to the number of potential partners in the network, (ii) closeness centrality, i.e., the shortest path from one node to all other nodes in the network, and (iii) betweenness centrality, i.e., the importance of a node for connecting different parts of the network. Centrality measures were computed using the software Pajek 4.10 (Batagelj & Mrvar, 1998). To test whether the node-level parameters varied significantly between the species-based network and the constrained trait-based network, we performed overall generalized linear models (GLM) and linear regressions (according to the distribution of the metric values).

Interacting partners and phenotypic traits

We recorded a total of 272 individual bee-flower interactions within our community of 10 plant and 28 bee species (Tables S1, S2). For plant traits, nectar holder depth (range = 0.98–11.63 mm) of sampled flowers (n = 131 flowers from the six out of 10 plant species with measurable nectar holder depths) was positively correlated (mean R² = 0.516 ± 0.101) with maximum floral display size (range = 7.0–53.44 mm) for most (four out of six) of the plant species subsampled for trait measurement (Fig. S2A). For bee traits, intertegular distance (range = 1.28–5.47 mm) was positively correlated (R² = 0.817) with proboscis length (range = 1.19–7.86 mm) across the subset of individuals (this included those collected during this experiment and some additional bees from a private collection held by T. Ings) measured in the laboratory (Fig. S2B).

Species-based network versus constrained functional trait-based network

The species-based network (Fig. 1A) comprised 10 plant species from four families (Table S1) and 28 bee species from six families (Table S2). Predicted NHD of flowers varied from 1.00 (open flowers) to 10.00 mm, whereas the intertegular distance of bees varied from 1.37 to 6.42 mm. In the nodes of the equivalent functional size-based network, where the number of nodes was constrained to match that of the species-based network (Fig. 1B), the number of species grouped within nodes ranged from 1 to 4 species in the case of the flowers visited (Table S3), and from 1 to 9 species in the case of the bees (Table S4). See also Fig. S3 for the distribution of the interacting flowers and bees in their respective size-based clusters. The interaction matrices (Fig. S4) and network representations (Fig. 1) of the two networks depict a regular pattern of flower-bee interactions that reveals the important role of body size. This pattern is more evident in the functional size-based network, where all intraspecific variation in body size has been considered, and size ranges do not overlap among nodes. In this network, a general pattern of size matching can be observed between the NHD of the flowers and the ITD of the bees (see Fig. 1B, where nodes have been ordered by size).

The two networks were quite similar in their structural parameters, although the species-based network showed lower values in C_q, IE and WNODF (Table 1 and Table S5). However, the species-based network was more specialized (H₂′) compared to its functional size-based counterpart (Table 1). By changing from the taxonomic to the functional trait-based representations, nestedness was the most affected parameter, jumping from 17.38 in the species-based network to 33.32 (Table S5) in the functional size-based network. The two webs showed significantly lower values of connectance (C_q) and interaction evenness (IE) than the randomly assembled networks (Table 1). However, both networks were significantly more specialised (H₂′) than the networks generated by the null model (Table 1). At the node level, bees were significantly more specialized in the species-based network compared to the functional size-based network (F_1,54 = 13.43, P < 0.001). There was no significant variation in the specialization degree of plants between the networks, and strength did not vary among networks for either bees or plants (Table 1).

Centrality varied between the two networks, as shown by the linear model fitted for the closeness metric in the case of bee-nodes (F_1,54 = 22.39, P <0.001). The mean value of closeness in the species-based network (0.36 ± 0.06) was significantly lower than that of the functional trait-based network (0.42 ± 0.07). While B. lucorum stands out for its central position in the species-based network (Table S6 : normalized degree = 0.60, closeness = 0.47, betweenness = 0.14), the node with highest value of centrality was B13 (normalized degree = 0.60, closeness = 0.49, betweenness = 0.12) for the size-based network. B13 groups individuals of three different species of Apidae (A. mellifera, B. lapidarius and B. pascuorum) with intertegular distances ranging from 3.69–3.83 mm (Table S4). For the plants, Centaurea nigra showed the most central position (normalized degree = 0.46, closeness = 0.46, betweenness = 0.31) in the species-based network, whereas the node F03 was the most central in the functional trait-based network (normalized degree = 0.82, closeness = 0.64, betweenness = 0.52). This flower node (F03), with a NHD ranging from 3.21 to 4.15 mm (Table S3), comprised flowers of two plant species: C. nigra and Leucanthemum vulgare.

Unconstrained functional trait-based network

When the functional size-based network was constructed with the optimal number of clusters according to trait variation (Fig. 2), it comprised five flower nodes that encompassed one to four plant species (Table S7), and four bee nodes that covered seven to 17 species (Table S8). Figure S5 shows the distribution of the interacting flowers and bees in their respective size-based clusters. A general pattern of size matching between bees and flowers was also evident in this unconstrained functional size-based network (Fig. 2). Thus, bees with longer proboscises tended to interact more frequently with flowers with deeper nectar holder tubes, whereas bees with shorter proboscises tended to interact more frequently with flowers with shallow nectar holder tubes.

Like the constrained functional size-based networks (with equal number of nodes as interacting species), the unconstrained network was significantly more specialized than the randomly assembled networks of the null models (Table 1). It also had significantly lower values of C_q and IE than the randomly assembled networks, a pattern also observed for the constrained network. The observed value of nestedness did not differ significantly from the null models (Table 1, Table S5).

At the node level, the centrality metrics (see Table S9) show that among the plants the node F02, composed of C. nigra and L. vulgare, with a NHD ranging from 3.02 to 4.15 mm (Table S7), occupied the most central position. Among bees, those with highest values of centrality were B02 and B03 (both with normalized degree = 1.00, closeness = 0.73, betweenness = 0.23). These nodes include bees belonging to four different families that show intermediate sizes (Andrenidae, Apidae, Megachilidae and Melittidae), and cover intertegular distances ranging from 3.20 to 3.83 mm and from 3.88 to 5.05, respectively (Table S8).