Individual-based plant–pollinator networks are structured by phenotypic and microsite plant traits

2832 | Journal of Ecology. 2021;109:2832–2844.wileyonlinelibrary.com/journal/jec Received: 6 November 2020 | Accepted: 4 May 2021 DOI: 10.1111/1365-2745.13694 RESEARCH ARTICLE Individualbased plant– pollinator networks are structured by phenotypic and microsite plant traits Blanca ArroyoCorrea1 | Ignasi Bartomeus1 | Pedro Jordano1,2 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. Journal of Ecology published by John Wiley & Sons Ltd on behalf of British Ecological Society 1Departamento de Ecología Integrativa, Estación Biológica de Doñana, EBDCSIC, Seville, Spain 2Departamento de Biología Vegetal y Ecología, Universidad de Sevilla, Seville, Spain Correspondence Blanca ArroyoCorrea Email: [email protected] Funding information Consejo Superior de Investigaciones Científicas, Grant/Award Number: JAEINT18_EX_0080; Ministerio de Ciencia e Innovación, Grant/Award Number: CGL2017082847P, CGL201792436EXP and FPU19_02552 Handling Editor: Ryan Phillips Abstract 1. The biotic and abiotic context of individual plants within animalpollinated plant populations can influence pollinator foraging behaviour. Pollinator movements regulate pollen flow among plant individuals, and ultimately determine individual plant reproductive success. Yet the underlying drivers of this context dependency of interactions at the population level and their functional consequences for individuals remain poorly known. 2. Here we used a wellcharacterised population of Halimium halimifolium (Cistaceae), a Mediterranean shrub species, in combination with exponential random graph models (ERGMs) to evaluate how the intrapopulation variation in plant attributes configures individualbased plant– pollinator networks and determines their reproductive outcomes. Specifically, we assessed (a) how the intrinsic (i.e. phenotype and phenology) and extrinsic (i.e. microsite) plant attributes influenced the emerging configuration of the bipartite plant– pollinator network and the unipartite plant– plant network derived from pollinator sharing, and (b) how these plant attributes combined with the network topological position of individual plants affect their female fitness, measured as the total seed weight per plant. 3. We found that both intrinsic and extrinsic plant attributes contributed substantially to explain the configuration of both the bipartite and the unipartite pollination network. Besides the effects of plant attributes, the functional group to which pollinator species belonged was also important to determine the variance in plant– pollinator interaction odds, while the probability of plants to share more pollinator species was additionally influenced by the spatial distance between those plants. Furthermore, our results showed that these influences of plant attributes on network structure can be translated into functional outcomes at the plant individual level, with direct consequences for intrapopulation fitness variation. 4. Synthesis. This study builds towards a better understanding of the multiple drivers underlying the context dependency of plant– pollinator interactions and how they mediate the reproductive outputs of individual plants within a population. The application of our analytical framework allows a conceptual shift from descriptive to | 2833 Journal of Ecolog y ARROYOCORREA Et Al. 1 | INTRODUCTION The foraging movements of pollinators influence how pollen transfer occurs among animalpollinated plant individuals within a population (Morris et al., 1995), ultimately influencing their reproductive outputs. Within a generalist pollination system, a plant population is usually composed of individuals differing in their level of pollinator attraction, and therefore, the pollinators visiting a given plant are a subset of the pool of potential pollinators (Gómez et al., 2007; Herrera, 2005). These interindividual differences in the pollinator assemblage of individual plants partially contribute to generate nonrandom patterns of pollen transfer events (Valverde, 2017), which may affect population dynamics and the local structuring of genetic diversity through a network of assortative mating events. Individual plants within most populations differ in their phenotype and this variation among plants (e.g. in plant and flower size or flowering phenology) may cause differences in their pollination niche due to distinctive preferences by different pollinator guilds (Gómez et al., 2014). Therefore, individual plants displaying similar phenotypes are expected to interact with a similar pollinator assemblage, promoting mating events among those plants (Gómez et al., 2011). Furthermore, dissimilarities in phenological patterns among pollinator species may create temporal variation in pollinator availability, leading to differences in the pollinator assemblage interacting with plants diverging in flowering phenology (Fox, 2003; Primack, 1985). On one hand, individual plants with higher flowering synchrony within their population are expected to share pollinators with more conspecifics; but, on the other hand, they may also interact with fewer pollinators due to intraspecific competition. Thus, a higher flowering synchrony among conspecifics can influence the level of pollinator sharing and therefore, increase intrapopulation assortativity in mating (Weis, 2005). Variation in the pollinator niche among individual plants may also be influenced by their spatial position in the population (RodríguezRodríguez et al., 2017). The spatial position defines the ecological context of plants, such as the microtopographical conditions and the local plant community to which they are exposed (e.g. intraspecific and interspecific competition), generating different aspects of context dependency of interaction outcomes (Chamberlain et al., 2014). Because pollinator species may tend to forage in specific microsites, the finescale location of individual plants within a population is expected to affect their interactions with pollinators and consequently, the patterns of pollinator sharing among conspecific plants (RodríguezRodríguez et al., 2015). Most plant– pollinator interaction studies do not address this spatial component, which might entail important constraints given the finescale division of resources among insect pollinators in a plant population (Dupont et al., 2014) and the spatial heterogeneity in neighbourhood composition (Janovský et al., 2013). Plant– pollinator communities have been widely studied during the last decades using bipartite interaction networks (Bascompte & Jordano, 2007), where pollinator species are connected to the plant species they visit. However, speciesbased networks overlook the existing intraspecific variation by pooling all the individuals in a single group (Dupont et al., 2014; Tur et al., 2015). In contrast to the communitylevel approach, individualbased networks can be used to unveil processes at the population level (e.g. RodríguezRodríguez et al., 2017; Soares et al., 2021; Valverde et al., 2016) and convey a direct link to analyse demographic and evolutionary questions. In a bipartite, individualbased network within a plant population, pollinator species are connected to the individual plants they visit. When it is projected into a unipartite network, individual plants are visualised as nodes and the links connecting those nodes depict the level of pollinator sharing, which may serve as a proxy of mating probabilities (Fortuna et al., 2008). Thus, the unipartite network generated with this procedure is undirected, in contrast with real pollentransfer networks. This undirected unipartite network allows us to elucidate how interactions and hence, potential mating events, are structured within a plant population. Because intrapopulation variation in phenotype, phenology and microsite may determine how individual plants interact with pollinators, it is expected to also affect the relative position of those plants within the unipartite network. While there are clear hypotheses suggesting that plants occupying central positions may benefit from a higher level of pollinator sharing with conspecifics, increasing the probability of mating events and hence their reproductive success (Gómez & Perfectti, 2012; Soares et al., 2021), it is not trivial how to translate network topologies into functional implications. Recently developed exponential random graph models (ERGMs) allow us to evaluate how plant attributes influence the configuration of complex networks, advancing their study from descriptive metrics into a more cohesive predictive framework (Harris, 2014). Such analyses are useful to test hypotheses about processes that shape networks by modelling the probability of interaction establishment as a function of the characteristics of individuals (nodes) and the value of other links within the network (Snijders et al., 2010). ERGMs differ from traditionally used descriptive approaches in that it encourages and supports identifying underlying mechanisms that explain the whole structure of the network. Although this analytical tool has been widely used in social network studies in recent years, it predictive research on the evolutionary and ecological processes that give rise to complex ecological networks at the population level. KEYWORDS Doñana, ecological networks, Halimium halimifolium, intraspecific variation, mutualism, plant fitness, plant– pollinator interactions 2834 | Journal of Ecology ARROYOCORREA Et Al. has rarely been considered for ecological network analysis (but see Miguel et al., 2018). Here we aim to investigate how the intrapopulation variation in intrinsic (i.e. phenotype and phenology) and extrinsic (i.e. microsite) plant attributes shapes the structure of individual plant– pollinator networks and influences their functional outcomes (Figure 1). To address this question, we recorded individual plant– pollinator interactions in a population of the Mediterranean shrub Halimium halimifolium (Cistaceae). First, we assessed how the intrinsic and extrinsic plant attributes influenced the emerging configuration of the bipartite plant– pollinator network and the unipartite plant– plant network derived from pollinator sharing. Second, to evaluate how this association between plant attributes and network structure translates into functional outcomes, we analysed how individual female plant fitness (i.e. the total seed weight per plant) was affected by the combination of plant attributes and the topological position of individual plants within the pollination network. 2 | MATERIALS AND METHODS 2.1 | Study site The study was performed in Doñana National Park (37°07′52.2″N 6°31′40.9″W, 12 m a.s.l.) within a 1.2 ha area bounded on its four sides by a stream, pine forests and seminatural grasslands. Doñana National Park is located on the Atlantic coast of southwestern Spain, in an area with a Mediterraneantype climate where the vegetation is composed mainly of sclerophyllous shrublands. Halimium halimifolium (Cistaceae) is a Mediterranean shrub that occupies the slopes of stabilised sand dunes (DíazBarradas & GarcíaNovo, 1987) and usually blooms in late spring (May– June; Herrera, 1988) when all other cooccurring plant species have finished their flowering season. The study site is largely dominated by this plant species, which represents nearly 60% of the total vegetation cover. Although other plant species are also present in our study area (Figure S2), H. halimifolium was the only species in flower during the study period. Hence, we did not consider competition or facilitation for pollinators among the different plant species. Each H. halimifolium plant can produce up to 1,000 hermaphrodite yellow flowers (Figure S1) at once and flower opening occurs synchronously each day in the population. Flowers usually exhibit a conspicuous dark brown spot at the base of each petal, being the presence and size of these spots highly variable among individuals within populations (Figure S1). Although these dark spots have not yet been proven to act as floral guides in H. halimifolium, this kind of petal marks is widely recognised to enhance pollinator attraction in multiple plant species from a wide range of families (de Jager et al., 2017). As this plant species is selfincompatible and cannot reproduce asexually, crosspollination promoted by insect visitation is necessary to complete seed production and allow population persistence (Talavera et al., 1997). 2.2 | Sampling We conducted surveys to record pollinator visitations in the study site during the peak flowering period of H. halimifolium (25 days between May and June 2019). Prior to data collection, all individual FIGURE 1 Conceptual diagram illustrating the relationship between the bipartite and the unipartite individual plant– pollinator networks (a) and its consequences for plant reproductive success (b). We use a minor, simplified, set of plant traits to illustrate that individual plants differ in phenotypic traits (e.g. size, number of flowers and presence of flower guides) and microsite characteristics (e.g. different ground colour tones). (a) Links in the bipartite network represent the number of visits from each pollinator species to each individual plant (top), while links in the unipartite network depict the number of pollinator visits shared between each individual plant pair (bottom). (b) The number of offspring per individual plant is expected to increase with the number of pollinator visits shared with other individual plants, as it increases the chances of being involved in pollen transfer and, eventually, in mating events. Pollinator icons were obtained from www.divul gare.net [Colour figure can be viewed at wileyonlinelibrary.com] (a) (b) ABCD Number of pollinators shared with conspecifics Number of offspring per plant A C D B | 2835 Journal of Ecolog y ARROYOCORREA Et Al. plants of the population (160 individuals) were labelled. Interactions with pollinator species were recorded in 60 of those individuals, which were selected using a stratified random sampling where both isolated and aggregated plants were well represented (Figure S3). Along the study period, we performed four 30min surveys on each individual plant using video cameras (GoPro HERO7, GoPro Inc.). Therefore, each of the 60 labelled plants was surveyed for 120 min across the flowering season, and thus, the entire survey lasted 7,200 min. Prior to each of these four surveys per individual plant, we recorded the total number of open flowers. To avoid withinday temporal bias in the pollinator assemblage sampling, these four surveys per plant were carried out at different times during the time period with open flowers (8:30 to 13:30), spanned out across four different sampling days. This sampling period encompasses the time of the day where the maximum activity of flower visitors to H. halimifolium occurs (Talavera et al., 1997). We assessed the completeness of sampling effort by calculating Chao asymptotic estimators (see Appendix S1 and Figure S4). Video cameras were placed approximately 30 cm away from the plant in each survey and therefore the number of flowers in view [35.21 ± 5.40 (mean ± SE)] depended on the flower abundance of those plants. To capture the variability in pollinator visits associated with the heterogeneous distribution of flowers within plants, cameras were set up at varying heights in each of the four surveys (range 20 to 180 cm tall). Video recordings were processed to annotate the identity of pollinator morphospecies visiting all the flowers in view and the number of times they interacted with those flowers (i.e. frequency of interactions). Pollinators were considered as all those insects landing on the flower and touching its reproductive structures. To facilitate pollinators' identification, a single representative specimen of each pollinator morphospecies present in the study site was captured in the field when possible, identified and vouchered at Estación Biológica de Doñana. Out of all pollinator morphospecies, 12 were identified at the species level, seven at the genus level, and only two at the family level. Independent surveys in our study area show that only one to three species with very similar morphology and behaviour occur within each of the lumped genera or families (Table S3). Besides, for the pollinator genera that performed most of the interactions, that is, Lasioglossum spp. (Hymenoptera: Halictidae), Dasypoda spp. (Hymenoptera: Melittidae) and Mylabris spp. (Coleoptera: Meloidae), only one species was identified with these independent surveys. For simplicity, we refer to these pollinator morphospecies as species hereafter. In order to characterise the plant– pollinator interaction assemblage, we first constructed a weighted bipartite network linking each individual plant with the pollinator species visiting its flowers. This individualbased network was built by creating an adjacency matrix A, where elements aij represent the number of interactions between the pollinator species i and the individual plant j. The number of interactions was calculated as the total number of times a pollinator species visited any flower of a given plant, and hence interaction data were analysed as the total number of visits. Second, to depict the pattern of shared pollinator species among plants, we generated the unipartite projection of matrix A by calculating the number of interactions shared by any two individual plants. The unipartite projection of matrix A for the P plants contains the parallel minimum of the observed interactions between two individual plants, so that the level of pollinator sharing is driven by the lowest communality in interactions with pollinator species. In this unipartite network, two individual plants are linked if they share at least one interaction with the same pollinator species. Hence, individual plants are not connected by planttoplant movements of individual pollinators, but by shared pollinator species. We assumed that sharing pollinator species is a surrogate indicator of the potential for effective pollen transfer among plants, given that most pollinators belonging to the same species have similar foraging behaviour (Fenster et al., 2004; Gómez & Perfectti, 2012). Thus, pollinator individuals from the same species are likely to show the same visitation patterns with individual plants. Two plants that are strongly connected in the unipartite network, and thus share most pollinators, will be more likely to be involved in mutual pollen transfer, and a potential mating event, than two plants with very different pollinators, connected by a weaker link. We thus considered the link weights estimated in this way as a proxy for the probability of mating events mediated by a given pollinator species (see Gómez et al., 2011; RodríguezRodríguez et al., 2015). 2.3 | Plant attributes and fitness estimation We obtained a set of intrinsic and extrinsic attributes for each surveyed individual plant, but we exclusively retained in the analysis those which were not redundant after testing for multicollinearity using VIF factors (threshold = 2, vif r package, Lin, 2012; for information about the excluded plant attributes see Appendix S1). The extrinsic attributes selected included microsite traits related to microtopographical conditions and competition for abiotic factors, which largely influence both pollinator visitation and the availability of abiotic resources to complete fruit production: (a) cover (m2) of intraspecific and interspecific neighbours in a 1 and 2 m radius, (b) distance to the stream, (c) distance to the habitat edge and (d) distance to the nearest tree (m). Intrinsic attributes comprised phenotypic and phenological traits associated with the ability of plants to attract pollinators: (a) maximum plant height (m), (b) total number of flowers during the flowering peak, (c) corolla diameter, that is, flower size (mm), (d) flower guide size (mm) and (e) flowering synchrony calculated by modifying the index by Marquis (1988) (Appendix S1). The specific biological predictions for the effects of these intrinsic and extrinsic variables on the bipartite and unipartite networks and the female plant fitness can be found in Table S1. Extrinsic variables and plant size were estimated using aerial images, an orthomosaic and a 3D surface model obtained with the performance of drone flights in the study population (see Appendix S1 for detailed methods of spatial data acquisition). To quantify the female reproductive success per individual plant, we first calculated the fruit set as the proportion of flowers setting fruit. To do this, we collected 10 inflorescences per plant and 2836 | Journal of Ecology ARROYOCORREA Et Al. counted the total number of flower buds produced initially and the number of fruits at the end of the season. Second, we sampled five fully developed fruits and counted the number of seeds per fruit to obtain the average seed production per fruit in each plant, and also measured the average seed mass with a precision scale. An overall female reproductive success (‘fitness’ hereafter) per plant was estimated by weighing the total number of seeds per plant by seed mass (Appendix S1). 2.4 | Statistical analyses 2.4.1 | Ecological correlates of the bipartite plant– pollinator network structure To analyse the ecological variables producing the overall structure of the individualbased pollination networks we built ERGMs (Lusher et al., 2013). An ERGM is a true generative statistical model of network structure and characteristics (i.e. it allows us to describe how data are generated in terms of probability): Node traits and local structural properties can be used to predict properties of the entire network (e.g. diameter, degree distribution, etc.). The probability of a given network, which can be reexpressed in terms of the conditional logodds of a single link between two nodes, is modelled as a function of terms that represent network, node or link features. The ERGM design is analogous to a classical generalised linear model (GLM) and implements a Markov chain Monte Carlo maximum likelihood parameter estimation. The use of these models for weighted networks allows us to test hypotheses about the underlying mechanisms shaping networks by modelling how the number of links (here interactions between a pollinator species and an individual plant or the number of pollinator species shared among individual plants) is affected by specific predictor variables. These predictor variables can be nodebased, dyadic and local structural covariates. Nodebased covariates model how node attributes affect the number of links they establish (e.g. plant size). Dyadic covariates explain differences in link values between two nodes as outcomes of other kinds of relationships between those nodes (i.e. link attributes, e.g., spatial distance between plants). Structural covariates, also called ‘network statistics’, describe aspects of network topology that are expected to affect link formation, being the total number of link weights the most basic structural term in weighted ERGMs. Local link configurations (e.g. motifs, Milo et al., 2002) that are hypothesised to occur more often or less often than expected by chance can be also included as structural covariates. Yet we focused only on the total number of link weights as the unique local structural predictor. The structural term representing the total number of link weights mentioned above, referred to as ‘sum’ in ERGMs, models the general propensity of nodes to be connected to other nodes, and it is equivalent to having an intercept within a GLM. Therefore, this framework allows us to account for the fact that link values can be dependent on the values of neighbouring links or some other aspect of network topology, directly dealing with the nonindependence related to this. The simplest ERGM controls for endogenous effects originated only from the pattern of links among nodes under a given reference distribution and hence, this model only incorporates structural covariates (Morris et al., 2008). However, the number of links between two nodes is expected to depend not only on the distribution of links joining two nodes, but also on exogenous effects of node and link attributes. A key choice in estimating weighted ERGMs is the selection of a reference distribution, which specifies a process underlying the count data before ERGM predictor terms are added. Because our plant– pollinator networks are based on counts of interactions, we used a Poisson structure, which effectively captures the underlying link distribution. The underlying Poisson distribution of the weighted ERGM allows the interpretation of the coefficients as the increase or decrease in the number of times that each plant– pollinator pair interacts (logged). We first fitted a Poissonreference ERGM with the structure of the bipartite plant– pollinator network as the response variable, the sum of the link weights as a structural covariate (i.e. endogenous effect), and the intrinsic and extrinsic attributes of individual plants as nodebased covariates (i.e. exogenous effects). The power of ERGMs is that we can also incorporate into the model pollinator nodebased attributes that can influence the odds of establishing a link, in this case, the pollinator functional group (see Table S2 for pollinator classification into functional groups). In this ERGM, the expected number of interactions between pollinator species and individual plants in the network is modelled as a random variable with a Poisson distribution. By associating the structure of the bipartite network to specific plant attributes we aimed to assess the contribution of each predictor variable of plant and pollinator attributes to the overall plant– pollinator network configuration. 2.4.2 | Ecological correlates of the unipartite plant– plant network structure To assess the association between plant attributes and the unipartite network structure depicting the overall patterns of pollinator sharing among conspecifics, we performed a second Poissonreference ERGM. We modelled the structure of the unipartite plant– plant network as a function of the sum of the link weights as a structural covariate, and the dissimilarity in intrinsic and extrinsic plant attributes between individuals as nodebased covariates. For each plant attribute, the dissimilarity is explicitly measured by the model as the absolute difference in two plants' values for that attribute. By doing so, we aimed to investigate whether individual plants with similar plant attributes shared more pollinator species, potentially leading to assortative mating events. We would expect that the higher the trait similarity among two plants, the greater their probability to share pollinator species. As a consequence of optimal foraging strategies of pollinators (Pasquaretta et al., 2017), they tend to move at short distances and thus, closer plants are more likely to share more pollinators than further plants and generate assortative matings. To account for this spatial effect, in this ERGM we also included | 2837 Journal of Ecolog y ARROYOCORREA Et Al. a matrix with the spatial distance between individual plants as a dyadic covariate. After performing both ERGMs, we compared the relative contribution of extrinsic and intrinsic plant attributes to explain link weights (i.e. number of interactions established) by expressing each predictor estimate as a percentage of the total sum of the absolute values of the estimates. All predictor variables were standardised to allow meaningful comparisons among them in relative importance. ERGMs were computed using sand and ergm r packages (Handcock et al., 2018; Kolaczyk & Csárdi, 2017). 2.4.3 | Functional consequences After assessing how extrinsic and intrinsic plant attributes shape the structure of the plant– pollinator networks, we evaluated how these influences translate into functional outcomes. The effects of plant attributes on plant fitness were tested using a GLM with a quasibinomial error distribution. We included the normalised plant fitness, which indicates how the fitness of a given individual departs from the population mean, as the response variable and the same extrinsic and intrinsic plant attributes as above as predictors. To account for the topological position of individual plants in the network, we estimated their centrality (Freeman, 1978) as a surrogate of the topological position of each individual. Centrality measures how well connected a given individual is with the rest of the cooccurring conspecifics owing to sharing pollinator species, and therefore its potential role as a pollen source and/or receiver. We assumed that a higher fraction of shared pollinator species reveals a higher potential for pollen transfer among two plants relative to other plant pairs with different pollinators. We estimated closeness centrality from the unipartite weighted network using the bipartite r package (Dormann et al., 2008). Closeness centrality is positively related to the shortest number of direct and indirect interactions between one node and all other nodes in the network (shortest paths; Freeman, 1978). In the context of this study, closeness is based on the reachability of an individual plant within the network, indicating how well pollen originating from a focal plant can reach any other plant or how well conspecifics' pollen can reach this focal plant. By including both plant attributes and closeness centrality as predictor variables in the GLM, we accounted for both the networkmediated effects and the external influences of plant attributes on fitness. Prior to model fitting, we tested for spatial autocorrelation of plant fitness and intrinsic plant attributes with Moran's I (spdep r package; Bivand & Wong, 2018). The relative contribution of both intrinsic and extrinsic plant attributes and plant centrality to explaining variation in fitness was assessed using the relaimpo r package (Grömping, 2006). This package uses the fitted linear model explained above with fitness as the response variable and intrinsic and extrinsic plant attributes as predictor variables and estimates the R2 contribution of those variables to explaining variation in fitness averaged over resampled orderings among regressors. All analyses were performed using R software version 3.5.3 (R Core Team, 2019). 3 | RESULTS We counted a total of 33,262 flowers from all surveyed individual plants [554.37 ± 108.64 (mean ± SE) per individual plant] across the study period. We recorded 2,774 interactions between 60 individual H. halimifolium plants and 21 pollinator species (Table S2). Most interactions were performed by Lasioglossum spp. (42.11%), Dasypoda spp. (18.78%), Apis mellifera (13.27%) (Hymenoptera) and Mylabris spp. (11.46%) (Coleoptera). The level of pollinator sharing (i.e. the number of visits by shared pollinator species) of each individual plant with conspecifics was 530.00 ± 56.73 (mean ± SE), while each individual plant shared pollinator species with 47.23 ± 1.36 conspecifics in the population. We observed substantial variation among individual plants in both extrinsic and intrinsic attributes (Table S4). The number of unique pairwise links between individual plants and pollinator species in the bipartite network was L(y) = 208 (Figure 2a), while the number of unique pairwise links among individual plants sharing pollinator species in the unipartite projection of the network was L(y) = 1,478 (Figure 2b). 3.1 | Ecological correlates of the bipartite plant– pollinator network structure We found a ‘sum’ effect in the ERGM fitted with the bipartite network as the response variable (Table 1), indicating that the probability of observing this network configuration was partially influenced by the total number of interactions among individual plants and pollinator species. However, additional variables (i.e. exogenous), such as plant attributes, were needed to explain the observed structure of the network. The model coefficients for exogenous variables (Table 1) represent the effects of plant or pollinator attributes on the likelihood of any plant– pollinator pair to increase the number of interactions. They indicate the expected change in the logged number of interactions between any plant– pollinator pair increasing (or decreasing) as a function of a given oneunit change in a specific plant or pollinator attribute. For example, our model showed that increasing the number of flowers in plants increased the odds of interacting with more pollinator species by a factor of exp(0.259) = 1.295, or nearly 30% (Table 1). Flower guide size increased the odds of plants to interact with more pollinator species, while these odds decreased with plant height and flowering synchrony. The direction of the effects of both intraspecific and interspecific neighbourhood cover on plant– pollinator interaction probability depended on the spatial scale considered (1 vs. 2 m radius). The effect of the cover of intraspecific neighbours was positive within a 1 m radius and negative within a 2 m radius. Meanwhile, the effect of the cover of interspecific neighbours was negative within 2838 | Journal of Ecology ARROYOCORREA Et Al. a 1 m radius and positive within a 2 m radius (see Table 1 for specific effect sizes). The functional group to which pollinator species belonged also had a significant effect on link establishment. The odds of establishing more plant– pollinator interactions especially increased when pollinator species were mediumand smallsized bees, by a factor of exp(2.114) = 8.281 and exp(3.480) = 32.460 respectively. Overall, intrinsic (i.e. phenotype and phenology) and extrinsic (i.e. microsite) plant attributes accounted for 9.54% and 8.01%, respectively, of the variation in plant– pollinator interaction odds, while the pollinator functional group explained 74.42% of this variation. 3.2 | Ecological correlates of the unipartite plant– plant network structure The ‘sum’ effect was significant in the ERGM with the unipartite network as the response variable (Table 2). Thus, the number of interactions between a given pair of plants (i.e. number of pollinator visits shared) partially depended on the total number of interactions in the unipartite network. The odds of two plants to share more pollinator species were positively associated with the dissimilarity between those plants in flower size, total number of flowers, distance to stream and cover of interspecific neighbours in a radius of 2 m. These odds decreased with the dissimilarity in plant height, flower FIGURE 2 (a) Weighted bipartite network representing interactions between individual plants (green nodes) and pollinator species (orange nodes). The links between nodes depict flower visitation interactions while the width of the links refers to the strength of the interaction (i.e. number of interactions recorded). The layout of the network representation was created using an energyminimisation algorithm. Pollinator species represent: And.sp - Andrena spp., Ant.dim - Anthaxia dimidiata, Ant.sp - Anthrenus spp., Api.mel - Apis mellifera, Bom.ter - Bombus terrestris, Bom.maj - Bombylius major, Cet.aur - Cetonia aurata, Cha.ill - Chasmatopterus illigeri, Das.sp - Dasypoda spp., Eri.ten - Eristalis tenax, Hel. ruf - Heliotaurus ruficollis, Lag.ser - Lagorina sericea, Las.sp - Lasioglossum spp., Lob.aen - Lobonyx aeneus, Mal.sp - Malthodes spp., Myl.sp - Mylabris spp., Nit - Nitidulidae, Oed - Oedemeridae, Pht.sp - Phthiria spp., Xyl.can - Xylocopa cantabrita, Zon.fla - Zonitis flava. (b) Weighted unipartite network illustrating the pattern of shared pollinator species among Halimium halimifolium individual plants. Two individual plants are linked if they shared at least one pollinator species visiting their flowers. Link width is proportional to the number of interactions shared between two individual plants. Node size refers to the female plant fitness while node position indicates the spatial location (X– Y coordinates) of individual plants [Colour figure can be viewed at wileyonlinelibrary.com] 120 127 129 136 137 143 150 152 156 158 161 163 164 1641 1645 1646 1651 1652 1654 1656a 1656b 1657 1659 166 1660 1662 1663 1664 1665 1671 1672 1674 1675 1678 168 1680 1685 1686 1689 169 1692 1693 1694 1697 170 171 177 187 189 195 196 197 198 199 2 200 201 202 3 out And.sp Ant.dim Ant.sp Api.mel Bom.ter Bom.maj Cet.aur Cha.ill Das.sp Eri.ten Hel.ruf Lag.ser Las.sp Lob.aen Mal.sp Myl.sp Nit Oed Pht.sp Xyl.can Zon.fla (a) (b) | 2839 Journal of Ecolog y ARROYOCORREA Et Al. guide size, flowering synchrony, distance to the nearest tree and cover of intraand interspecific neighbours in an r = 1 m (Table 2). We also found that the odds of pollinator sharing between plants decreased with the spatial distance between those plants by a factor of exp(1.790) = 5.989. Intrinsic and extrinsic plant attributes explained 15.65 and 15.88%, respectively, of the variation in the odds of pollinator sharing, while the spatial distance between individual plants accounted for 25.58% of this variation. 3.3 | Functional consequences Plant fitness and intrinsic plant attributes were not spatially autocorrelated (Table S5). The GLM with individual plant fitness as the response variable revealed that the topological position of individual plants within the unipartite network (i.e. closeness centrality) was positively associated with plant fitness (Table 3; Figure 3a). We also found that individual plant fitness was associated with intrinsic and extrinsic plant attributes, after controlling for the networkmediated effects (Table 3). Interestingly, individual plant fitness increased with the total number of flowers and flowering synchrony (Figure 3b,c); TABLE 1 Summary of the exponential random graph model (ERGM) showing the effects of individual plant attributes and pollinator functional groups on the structure of the bipartite plant– pollinator network. The ‘sum’ effect assesses the influence of the plant– pollinator link weights' distribution on the configuration of the bipartite network. Estimates of predictor variables with significant effects indicate the expected change, increasing or decreasing, in the logged number of plant– pollinator interactions given a oneunit change in that predictor variable. As the first pollinator functional group (beeflies) is used as a reference, the interpretation of the effects of the remaining groups is relative to beeflies. Significant values (p < 0.05) appear in boldface Estimate ± SE z value p Sum −0.963 ± 0.149 −6.472 <0.001 Intrinsic plant attributes Maximum height (m) −0.144 ± 0.033 −4.375 <0.001 Flower guide size (mm) 0.204 ± 0.028 7.394 <0.001 Flower size (mm) 0.021 ± 0.032 0.674 0.500 Total number of flowers 0.259 ± 0.023 11.297 <0.001 Flowering synchrony −0.517 ± 0.049 −10.486 <0.001 Extrinsic plant attributes Distance to stream (m) −0.073 ± 0.055 −1.325 0.185 Distance to tree (m) 0.006 ± 0.035 0.169 0.866 Distance to edge (m) 0.02 ± 0.045 0.455 0.649 Cover of intraspecific neighbours in r = 1 m (m2) 0.215 ± 0.051 4.202 <0.001 Cover of interspecific neighbours in r = 1 m (m2) −0.317 ± 0.062 −5.115 <0.001 Cover of intraspecific neighbours in r = 2 m (m2) −0.109 ± 0.048 −2.283 <0.001 Cover of interspecific neighbours in r = 2 m (m2) 0.221 ± 0.041 5.426 <0.001 Pollinator functional group Hoverflies 0.211 ± 0.222 0.952 0.341 Large bees 0.223 ± 0.179 1.247 0.212 Large beetles 1.482 ± 0.145 10.208 <0.001 Mediumsized bees 2.114 ± 0.135 15.625 <0.001 Mediumsized beetles −0.933 ± 0.17 −5.485 <0.001 Small bees 3.480 ± 0.137 25.364 <0.001 Small beetles −0.484 ± 0.196 −2.472 <0.001 TABLE 2 Summary of the exponential random graph model (ERGM) showing the effects of individual plant attributes on the structure of the unipartite network derived from pollinator species sharing among conspecific plants. The ‘sum’ effect assesses the influence of the plant– plant link weights' distribution on the configuration of the unipartite network. Estimates of predictor variables with significant effects indicate the expected change, increasing or decreasing, in the logged number of pollinator species shared given a oneunit change in that predictor variable. Significant values (p < 0.05) appear in boldface Estimate ± SE z value p Sum 3.005 ± 0.028 108.984 <0.001 Intrinsic plant attributes Maximum height (m) −0.119 ± 0.012 −10.132 <0.001 Flower guide size (mm) −0.222 ± 0.01 −21.219 <0.001 Flower size (mm) 0.138 ± 0.011 12.335 <0.001 Total number of flowers 0.134 ± 0.006 20.620 <0.001 Flowering synchrony −0.482 ± 0.015 −31.978 <0.001 Extrinsic plant attributes Distance to stream (m) 0.197 ± 0.037 5.265 <0.001 Distance to tree (m) −0.047 ± 0.013 −3.645 <0.001 Distance to edge (m) −0.025 ± 0.013 −1.855 0.064 Cover of intraspecific neighbours in r = 1 m (m2) −0.112 ± 0.014 −8.116 <0.001 Cover of interspecific neighbours in r = 1 m (m2) −0.471 ± 0.017 −28.342 <0.001 Cover of intraspecific neighbours in r = 2 m (m2) 0.013 ± 0.015 0.875 0.381 Cover of interspecific neighbours in r = 2 m (m2) 0.246 ± 0.01 24.184 <0.001 Spatial distance between plants (m) −1.790 ± 0.186 −9.552 <0.001 2840 | Journal of Ecology ARROYOCORREA Et Al. and decreased with distance to the stream and cover of intraspecific neighbours in r = 2 m (Table 3; Figure 3d,e). Intrinsic plant attributes accounted for 50.64% of the variation in individual plant fitness, and extrinsic attributes explained 34.59% of this variation. Meanwhile, the topological position of individual plants within the network accounted for 14.78% of this variation. 4 | DISCUSSION Our findings revealed how the individual variation in plant attributes existing in natural populations gives rise to the overall structural patterns of individualbased plant– pollinator networks and its functional implications. We found that intrinsic (i.e. phenotypic traits and flowering phenology) and extrinsic (i.e. microsite) plant attributes had comparable influences on structuring both, the plant– pollinator bipartite network and the plant– plant unipartite network. Furthermore, our results also showed that network structure, and specifically the position of individual plants within this network, together with the direct effects of plant attributes, had considerable consequences for individual plant fitness. Overall, these findings suggest that the ability of individual plants to participate in shared interactions, and therefore to increase their reproductive success, may ultimately be determined by the intraspecific variation in microsite characteristics, phenotype and phenological traits. 4.1 | Ecological correlates of the bipartite plant– pollinator network structure Variation in traits of both plant and animal partners is increasingly recognised as forces structuring ecological networks at the level of interactions among species (Bartomeus et al., 2016). Similarly, intraspecific trait variation can play an important role in shaping individualbased networks, yet this aspect remains virtually unexplored in the literature (Miguel et al., 2018). By using ERGMs we were able to further recognise specific plant attributes (individual phenotype, phenology and microsite) that shape the observed individual plant– pollinator network structure. We thus assessed the relative contribution of extrinsic and intrinsic plant attributes to explaining the odds of interactions to increase for individual plants. Our findings showed that intraspecific trait variation in floral and vegetative phenotypes among plants influenced the establishment of interactions with pollinator species, which is in accordance with previous work (Dupont et al., 2014; Gómez & Perfectti, 2012; Weber & Kolb, 2013). We also demonstrated that this influence upscales to affect the structure of the individualbased pollination network. Specifically, we found that the probability of interacting with more pollinator species increased with the size of the flower guides, which is widely known to positively affect pollinator attraction (Waser, 1983). Although plant size is expected to increase plant attractiveness for pollinators, we found a negative effect of plant height on the probability to interact with pollinator species. This result is probably due to the nonuniform distribution of flowers in taller plants, as they usually have the highest flower density at the base of the plant. Contrary to our expectations, the individual flowering synchrony influenced the structure of the network by decreasing the odds of individual plants to interact with pollinator species. This result may respond to the fact that a high abundance of flowers in the population can decrease the chance of a pollinator to visit a given individual plant, leading to the socalled ‘dilution effect’ (Holzschuh et al., 2011), especially in highdensity stands. However, we found that the flowering synchrony of plants also increased their fitness. Our findings suggest that higher floral availability in the population can result in high pollen loads carried by pollinators, which can ultimately increase individual plant fitness. Another nonexclusive plausible explanation is that plants with higher synchrony can share pollinators with a larger set of conspecific plants, which may TABLE 3 Results of the generalised linear model (GLM) evaluating the effects of the network topological position and individual plant attributes on female plant fitness. The fitted GLM accounted for 84.81% of the variation. Fitness was estimated as the total number of seeds per individual weighted by the seed mass and normalised (See Materials and Methods: Plant attributes and fitness estimation). Significant values (p < 0.05) appear in boldface Estimate ± SE t value p Intercept −2.597 ± 0.121 −19.610 <0.001 Topological position Closeness 0.767 ± 0.215 3.576 <0.001 Intrinsic plant attributes Maximum height (m) −0.232 ± 0.157 −1.478 0.142 Flower guide size (mm) 0.127 ± 0.116 1.092 0.281 Flower size (mm) −0.103 ± 0.131 −0.783 0.437 Total number of flowers 1.501 ± 0.162 9.271 <0.001 Flowering synchrony 0.666 ± 0.237 2.808 0.007 Extrinsic plant attributes Distance to stream (m) −0.577 ± 0.219 −2.639 0.011 Distance to tree (m) −0.252 ± 0.171 −1.474 0.147 Distance to edge (m) 0.077 ± 0.221 −0.349 0.729 Cover of intraspecific neighbours in r = 1 m (m2) 0.259 ± 0.207 1.2503 0.218 Cover of interspecific neighbours in r = 1 m (m2) −0.059 ± 0.206 −0.286 0.776 Cover of intraspecific neighbours in r = 2 m (m2) −0.469 ± 0.192 −2.447 0.018 Cover of interspecific neighbours in r = 2 m (m2) −0.099 ± 0.202 −0.494 0.623