Drivers of Individual-based, Antagonistic Interaction Networks During Plant Range Expansion

2190 | Journal of Ecology. 2022;110:2190–2204.wileyonlinelibrary.com/journal/jec Received: 17 January 2022 | Accepted: 18 May 2022 DOI: 10.1111/1365-2745.13942 RESEARCH ARTICLE Drivers of individualbased, antagonistic interaction networks during plant range expansion Jorge Isla1 | Miguel E. Jácome- Flores1,2 | Daniel Pareja3 | Pedro Jordano1,3 This is an open access article under the terms of the Creative Commons Attribution-NonCommercial License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited and is not used for commercial purposes. © 2022 The Authors. Journal of Ecology published by John Wiley & Sons Ltd on behalf of British Ecological Society. 1Estación Biológica de Doñana, CSIC, Sevilla, Spain 2CONACYT- Centro del Cambio Global y la Sustentabilidad, Villahermosa, Tabasco, Mexico 3Departamento de Biología Vegetal y Ecología, Facultad de Biología, Universidad de Sevilla, Sevilla, Spain Correspondence Jorge Isla Email: [email protected] Funding information Cátedras- CONACYT, Grant/Award Number: 73215; Ministerio de Ciencia e Innovación, Grant/Award Number: C G L 2 0 1 7 - 8 2 8 4 7 - P , P R E 2 0 1 8 - 0 8 5 9 1 6 , S U M H A L , L I F E W A T C H - 2 0 1 9 - 0 9 - C S I C - 1 3 and POPE 2014- 2020 Handling Editor: Shurong Zhou Abstract 1. Range expansion in plant populations, especially at the colonization front, can be either limited by disproportionately large effects of antagonistic interactions or facilitated by their release. How the strength of antagonistic interactions changes along successional gradients during range expansion is still poorly documented, especially when diverse assemblages of plant antagonists (rodents, invertebrates and birds) combine within interaction networks. 2. We study the changes in individualbased, predispersal seed– pulp predator networks along a colonization gradient in a rapidly expanding Juniperus phoenicea population in Doñana National Park (SW Spain). Additionally, we analysed the role of individual plant traits and neighbourhood attributes in network configuration by using Exponential Random Graph Models. 3. Seven seed– pulp consumer animal species varied significantly in their frequency of interaction and prevalence. While invertebrate species were wellestablished in old and intermediately mature stands, greenfinch Chloris chloris was dominant at the colonization front. Variable species roles and spread of interactions among individual plants generated changes in the configuration of interactions during plant expansion. 4. Individual plant traits strongly determined the topology of these networks, although with differences between stands. Increasing individual crop size and seeds per cone increased the interaction odds of individual plants, while seed viability showed the opposite effect. The network topology at the colonization front appeared less driven by individual traits, possibly because of the short interaction history of this recently established area. The disproportionately large effect of C. chloris in these recently established stands, potentially resulted in large seed losses during range expansion. 5. Synthesis. Turnover of antagonistic interactions, characterized the colonization front, resulting in more heterogeneous interaction strengths among individual plants. We found no evidence for a complete or sizeable antagonistic release of J. phoenicea at the colonization front promoting this rapid expansion. It becomes necessary to explore interactions with seed dispersers to understand | 2191 Journal of Ecology ISLA et al. 1 | INTRODUCTION Fastpaced global change is currently impacting a number of crucial ecological factors that trigger plant range shifts. Recently, robustly supported evidences including spatial mismatches among interacting species caused by climate change (Schweiger et al., 2008), natural recolonization of abandoned agricultural lands (Escribano- Avila et al., 2014) or natural expansion (García et al., 2014) have shown plant abilities to naturally regenerate, disperse and colonize new habitats over short time spans. In most cases, a clear involvement of a diverse array of plant– animal interactions has been documented shaping the plant movements. For instance, mutualistic interactions such as animalmediated seed dispersal are necessary to trigger rapid responses for plants to respond efficiently to global change drivers (González- Varo et al., 2021). However, the role that antagonistic interactions may play in limiting or facilitating plant range shifts in response to these factors is largely unknown. Plant populations in natural range expansion scenarios are characterized by a mature population that acts as a source of propagules heading towards recently established areas, known as the colonization front (Shigesada & Kawasaki, 2002). Plant stand features such as individual density, spatial cover or demographic and genetic structure vary in these range shift processes (Excoffier et al., 2009; Petit, 2011; Shigesada & Kawasaki, 2002). Population expansion often involves a complex network of interactions with animals that either facilitate expansion (e.g. seed dispersers, pollinators) or constrain it (e.g. herbivores, seed predators). The outcome and balance between positive (mutualistic) and negative (antagonistic) interactions with other effects (e.g. interspecific competition) may pervasively determine the dynamics of plant movements (Svenning et al., 2014). It is important to explicitly represent interspecific interactions in forecasts of dynamic range expansion when most interacting species show correlated spatiotemporal trends in their effects and the number of interacting species is low (Svenning et al., 2014). The presence of active expansion fronts necessarily associates with such spatial and geographical variation in interaction conditions (Travis, 1996). One of the possible hypotheses linking natural expansion processes and plant– animal interactions is the enemy release hypothesis (Keane & Crawley, 2002). This release from enemies facilitates the expansion through a reduction or elimination of antagonistic agents, repeatedly documented during expansions of alien species (Meijer et al., 2016). Antagonistic plant– animal interactions related to reproductive fitness play a key role during processes restructuring plant populations. For example, by increasing interactions intensity or releasing their pressure, these could constrain or allow early regeneration processes (Keane & Crawley, 2002; Svenning et al., 2014) or determine species richness (Janzen, 1970). These are likely processes during colonization from source stands, yet we still have a limited understanding regarding how biotic interactions and ecological factors result in observed range dynamics. Antagonistic interactions with animals and their effect on plants can deeply affect plant reproductive output (Sallabanks & Courtney, 1992; Strauss & Irwin, 2004). These interactions could determine seed source limitation: plants are simply unable to produce enough propagules to ‘fill’ available target microsites for recruitment (Münzbergová & Herben, 2005; Schupp, 2002). Antagonistic and mutualistic interactions in nature are embedded, building complex networks with species from very different taxa (Fontaine et al., 2011; Morris et al., 2014) and result in combined effects on plant fitness (Harper, 1977). Predispersal seed predation occurs in advance to other potentially limiting effects that act on seed dispersal, such as recruitment or regeneration limitation (Nathan & Muller- Landau, 2000; Schupp, 2002). In fleshy fruited species, predispersal predation interactions may involve direct seed destruction (Fuentes & Schupp, 1998; González- Varo, 2010), fruit infestation by insect pests and other pathogens (García, 1998) and fruit/seed damage by pulp peckers and seed predators (Simmons et al., 2018; Snow & Snow, 1988). Subsequently, predation on dispersed seeds is performed by major seed predators and nonlegitimate seed dispersers including insects, mammals and birds (Hulme & Benckman, 2002; Janzen, 1971). The magnitudes of predispersal seed predation are usually relatively low (Janzen, 1971; Kolb et al., 2007; Xu et al., 2015), although sometimes can compromise a large part of the available propagules (Crawley, 2000; Guido & Roques, 1996). How assemblages of predispersal seed– pulp predators are reshaped along plant regeneration gradients remains underexplored (Sallabanks & Courtney, 1992). Tools from ecological networks theory have recently proven most effective in assessing the complex patterns of plant– animal interactions (Bascompte & Jordano, 2014). However, network models are typically built on speciesaveraged estimators, ignoring variability among individuals in their interaction patterns (Dupont et al., 2014; Melián et al., 2014). By averaging, the importance of individual biotic and abiotic context in the establishment of interactions how antagonistic and mutualistic plant– animal interactions balance during range expansion. Our study highlights the importance of an individualbased approach in understanding how interactions are structured and driven in natural changing landscapes. KEYWORDS antagonistic interactions, dispersal, ecological networks, landscape ecology, plant range shifts, plant– herbivore interactions, reproductive ecology, seedpredation 13652745, 2022, 9, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2745.13942 by Universidad De Sevilla, Wiley Online Library on [17/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 2192 | Journal of Ecology ISLA et al. is neglected (Rodríguez- Rodríguez et al., 2017; Thompson, 1988; Valverde et al., 2016). This approach has pervasive consequences for the inferences and results interpretation. For example, aggregated data from many individuals subject to spatiotemporal variation are used to produce specieslevel averages, which marginalize away the relevant (processlevel) scale (Clark et al., 2011). It is important to note that the partners in interactions in nature are individuals, not the species to which they belong. This fine scale is the most appropriate when assessing factors structuring interaction networks between plants and animals (Dupont et al., 2011). An individualbased approach is a powerful tool to understand the role of predispersal antagonistic interactions during plant range expansion, and factors that structure them. Any interaction shaped by individual phenotypes (with heritable potential) that could modify individual plant fitness can be an important selective force for trait evolution (Strauss et al., 2005; Strauss & Irwin, 2004). At the intrapopulation level, both intrinsic and extrinsic individual plant attributes play a role in interactions with predispersal antagonists (Schupp et al., 2019). The main intrinsic traits are those related to individual fecundity, plant physical features and seed– fruit traits. Meanwhile, extrinsic plant attributes such as location, isolation and neighbourhood also play a role in the assembly of antagonists (Schupp et al., 2019). Novel analytical tools like Exponential Random Graph Models (ERGMs) allow us to unravel the determining factors of the complex structure of interaction networks (Morris et al., 2008; Saul & Filkov, 2007) and are extremely effective in addressing individualscale variation. Analogously to a generalized linear model, ERGMs are based on a response variable which is the structure of the network itself (link distribution among nodes), and predictor variables with information associated with each specific network node (i.e. plant attributes). Thus, it is possible to compare the structure of an observed network with that generated by models that include or exclude specific node information. These models have been recently introduced in ecological research (Arroyo- Correa et al., 2021; Miguel et al., 2018), providing ways to infer causes determining the distribution of the interactions among network nodes (Kolaczyk & Csárdi, 2014). The relative roles of intrinsic and extrinsic traits in interaction network structures have been evaluated for mutualistic networks (Arroyo- Correa et al., 2021; Miguel et al., 2018; also see Gómez et al., 2011; Valverde et al., 2016), although how these factors drive antagonistic network topologies at the individual level remains unexplored. In this study, we focused on the predispersal antagonistic assemblage of seed– pulp predators of Juniperus phoenicea subsp. turbinata. Juniperus phoenicea is considered as a foundation species (Whitham et al., 2006) and has undergone a rapid expansion since the protection of Doñana National Park five decades ago (García et al., 2014). This plant range expansion setting, involving a diverse assemblage of predispersal antagonists, is ideal for examining how interaction networks are restructured at the individual level across colonization stages, and the factors that determine these shifts in network structure. We expect a reconfiguration of interaction networks along the colonization gradient that would explain the rapid expansion of this population (e.g. through a release of antagonistic interactions in recently established stands). Furthermore, we hypothesize that the individualbased network topologies are being driven by individual traits and neighbourhood context. Specifically, the main goals in this study were: (i) Describe the individualbased network of J. phoenicea and its assemblage of predispersal seed– pulp predators. (ii) Examine how these interactions are reshaped along a natural gradient of colonization through basic network descriptors and the species and the spread of interactions across individual plants. (iii) Assess the role of plant traits and neighbourhood context as drivers of network topology. (iv) Evaluate the consistency of network topology drivers during plant range expansion. 2 | MATERIALS AND METHODS 2.1 | Study site and study species The study area is located on the southwestern edge of the Reserva Biológica de Doñana (RBD) within the limits of the Doñana National Park (Figure S1). Until its protection in 1969, this area was managed by local inhabitants for wood harvesting, livestock and hunting (Granados et al., 1988). Since then, several plant species, notably J. phoenicea, have significantly expanded their distribution in the area by natural regeneration (García et al., 2014). Locally, its expansion process has produced a natural colonization gradient from mature stands to early colonization fronts, where J. phoenicea enters areas dominated by low sclerophyllous shrubland. Juniperus phoenicea subsp. turbinata (Guss.) Nyman (Cupressaceae) is a gymnosperm shrub inhabiting coastal dunes and rocky habitats in the western Mediterranean and in the Macronesian archipelagos (Adams, 2011). Species of the genus Juniperus are considered as foundation species in semiarid Mediterranean ecosystems (Whitham et al., 2006) and the juniperdominated woodlands on stabilized dunes are listed as Habitat of Community Interest (DOCE, 1992). Juniperus phoenicea is an anemophilous species characterized by very productive mast seasons occurring among years with very limited cone production (Jordano, 1993; Roques et al., 1984). Its reproductive cycle spans 2 years, and adult trees produce transformed fleshy cones also called galbules or arcestides (hereafter, cones). Cones are brownred coloured when ripe (Figure S2) with an average of five seeds per cone, consumed and dispersed by several thrush species and mediumsized generalist mammals, such as Meles meles and Vulpes vulpes (Jordano, 1993; Perea et al., 2012). Several invertebrate consumers have been recorded feeding on the cones and seeds all along the J. phoenicea distribution area (Roques et al., 1984; Turgeon et al., 1994; Ribes & Askew, 2009). The effects of each specific interaction are variable (Figure 1; Table S1). For example, the wasp Megastigmus amicorum (Hymenoptera: Torymidae) feeds on the seed embryo, while the larvae of the beetle Nanodiscus transversus (Coleoptera: Nanophyidae) feeds on the cone pulp. The antagonistic mammal assemblage is mainly 13652745, 2022, 9, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2745.13942 by Universidad De Sevilla, Wiley Online Library on [17/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 2193 Journal of Ecology ISLA et al. composed by mice Apodemus sylvaticus and Mus spretus (Rodentia: Muridae) that destroy the seeds to feed on the energyrich embryo of juniper species, although rodents can consume the outer pulp (García et al., 2000; Santos & Tellería, 1994). In the study area, larger mammal species including lagomorphs, carnivores and ungulates, may also play a role as seed predators (also as seed dispersers) of J. phoenicea, damaging a variable percentage of seeds when feeding on cones (Perea et al., 2012). However, we exclude these species since their role as seed predators occurs in late stages of recruitment. Finally, the only reported seed predatory bird species is the greenfinch Chloris chloris (Passeriformes: Fringillidae; Santos & Tellería, 1994). 2.2 | Sampling design Our sampling design consisted of three 1- ha plots along a natural regeneration gradient (Figure S1). These plots have all junipers individually identified and georeferenced. The mature stand is located in a mature juniper forest named ‘Sabinar del Marqués’ (MAR), with dominance of junipers and a high density of individuals (10,520 individuals/ha). The intermediate stand, named ‘Sabinar del Ojillo’ (OJI), is a dense juniper stand (9010 individuals/ha) also composed of mediterranean scrubland. The colonization front stand, ‘Sabinar de Colonización’ (COL) is characterized by scattered juniper individuals of recent formation and with the lowest density (2030 individuals/ha). We established five subplots [mean (±1 SD) of 522 ± 115 m2 subplots area] within each of the three main 1 ha plots, distributed regularly throughout the main plot. We used the individuals located in each subplot to randomly select 35 individual plants per stand (N = 105 plants in total). Focal plants were characterized and monitored for individual interactions during two fruiting seasons: October 2018– May 2019 and October 2019– June 2020. 2.3 | Sampling interactions We monitored each of the 105 focal plants during the two consecutive study seasons, using different methods (see Suppl. Mat. for details). To collect data on rodent interactions, we combined two sampling methods: cameratrap survey and livetrapping. We conducted eight direct nighttime trapping campaigns and the phototrapping sampling effort involved a total of 201,600 h of recording, both methods evenly distributed among focal plants. For both methods we consider the presence of a species under a focal plant as a potential interaction. We quantified invertebrate interactions by dissecting 50 mature cones from each focal plant per season (number of analysed cones = 7906 cones). We identified antagonistic species and their markings based on Roques et al. (1984). When a parasitoid was recorded (Figure 1), we assigned the event as a detection of the host species, and indicative of one interaction of the plant with the corresponding host pulp– seed predator species (Table S1). Finally, C. chloris interactions were quantified by seedtraps placed under focal plants to collect the remnants of cones falling during their visits. Based on video records of foraging finches during full visits to the plants [mean (±1 SD) of 4.5 ± 2.6 cones consumed per visit], we FIGURE 1 The predispersive consumers of Juniperus phoenicea pulp and seeds found in the Doñana Biological Reserve. In addition, the recorded parasitoid species or morphospecies are included. The lines starting from the consumers go to the cone pulp or seed depending on their specific target. Modified from Ribes and Askew (2009). Chalcididae sp. 2 Nanodiscus transversus Megastigmus amicorum Trisetacus q uadrisectus Elachista sp. Mesophleps oxycedrella Entedon juniperi Aprostocetus galbulus Eupelmus acinellus Braconidae sp. Chalcididae sp. 1 ? Chloris chloris Apodemus sylvaticus Mus spretus 13652745, 2022, 9, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2745.13942 by Universidad De Sevilla, Wiley Online Library on [17/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 2194 | Journal of Ecology ISLA et al. estimated that every four damaged cones found in the seedtraps conservatively indicate one visit event by C. chloris in an individual plant. We carried out 35 checks of the seedtrap during the entire sampling, in which we counted the number of cones attacked, and emptied the trap. To handle the same interaction units for all data, we calculated the frequency of occurrence of each animal interaction for its survey. For example, how frequently an insect species was present in a cone relative to all the analysed cones, or how frequently we detected at least one C. chloris visit during weekly seedtraps checks. 2.4 | Plant characteristics We thoroughly sampled a set of intrinsic and extrinsic plant traits known to influence animal preferences across the 105 focal plants (Table S2). We considered plant traits that could drive interactions with predispersal seed– pulp predators hierarchically (Sallabanks, 1993): first, the individual neighbourhood context, then general plant traits and finally at a higher scale of detail, cone quality traits. To characterize the juniper neighbourhood density and productivity of each plant we used the georeferenced location of all the individuals in each stand and estimated by direct count the cone production in all the juniper individuals growing in a buffer area of 100 m2 surrounding each focal individual. For each focal individual, we recorded its height, the two maximum diameters of its canopy projection, canopy area and total cone crop size. Direct counts of cones were carried out by scanning the whole plant canopy area and counting the cones with a handcounter. We used the harvested cones (50 cones per plant) to measure cone traits. For each plant, we measured the average values of maximum length and diameter of the cones, total fresh mass, pulp mass, seed mass, oneseed mass, number of seeds per cone and seed viability was estimated by flotation procedure (Table S2; Figures S 4 – S 6 ). 2.5 | Individualbased, antagonistic interaction networks of J. phoenicea We constructed four weighted, individualbased, bipartite networks between J. phoenicea focal plants and their predator species, one for each of the studied stands and a general one with all the data combined. Notice that our networks are partially individual networks (i.e. bipartite networks with two modes, plants and animals), with one of their modes (plant nodes) being represented by individual plants, but the other mode (animal nodes) represented by the antagonistic animal species. However, we use a phytocentric approach where the individualbased scale is only applied to plants and their individual assemblages. We used the three perstand networks to evaluate changes in network configuration and node strength distribution during plant range expansion. We pooled data collected from the three stands (105 focal plants) into an overall network that was useful as a ‘big picture’ to detect general trends in the effect of cone and plantlevel traits in the network topology. In addition, perstand networks were also used to assess amongstand variation in the role of individual traits and neighbourhood attributes on the detected network topologies across the colonization gradient. To detect structural changes in interaction organization between the three study areas a pool of network metrics was computed. To evaluate how the weight of interactions was spread out or concentrated in a few plants we calculated species strength. This metric was also calculated for predator species, to assess their role in each network. This nodelevel metric measures the importance of each node of one mode in the network from the perspective of the other mode nodes (Bascompte et al., 2006). In addition, as a measure of how extensively the interactions with each animal species are distributed among plant individuals in each stand, we calculated the prevalence of predator species (% plants with specific interaction recorded). At the network level, we calculated metrics that would allow us to detect changes in the overall network topology (weighted connectance, weighted nestedness and modularity). To assess the antagonistic load in each stand we calculated Weighted Connectance (wC), the number of links in the whole network relative to the total possible number of links (Bersier et al., 2002). To determine the amount of hierarchical structure in a network, we computed Weighted Nestedness (wNODF). High wNODF strengthens a core of highly connected nodes which, coupled with increased connectance of peripheral nodes, provides a cohesive structure to the network (Bascompte et al., 2003). To detect densely connected groups of nodes (e.g. animal species that tend to interact with the same subset of individual plants), with sparse connections to nodes in other groups, we calculated Modularity (Olesen et al., 2007). Both networklevel and nodelevel metrics were computed with r package bipartite, and in the case of modularity we use the Beckett algorithm (Dormann et al., 2008; Dormann & Strauss, 2014). To address network similarities, we computed the Hamming distance between network graphs, which relies on the minimum number of link additions/deletions required to mutate one network into the other (Butts & Carley, 2005). In this context, higher Hamming distance between networks A- B than between A- C indicates that the topology of A is more similar to C than it is to B. Statistical significance of network metrics was calculated by comparison with a null model, contrasting the observed networks with 100 randomly generated networks using null models which maintain network connectance (‘vaznull’ function in the r package bipartite, Dormann et al., 2009). We used bootstrapping with replacement to estimate 95% confidence intervals (CIs) for network metrics and to assess possible overlap in estimated metrics values between the three study networks. 2.6 | Assessing drivers of network topology To analyse the effect of individual and neighbourhood plant traits on network topology we estimated an Exponential Random Graph Model (ERGM; Kolaczyk & Csárdi, 2014; Lusher et al., 2013). Before fitting ERGM models on network structure and performing 13652745, 2022, 9, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2745.13942 by Universidad De Sevilla, Wiley Online Library on [17/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 2195 Journal of Ecology ISLA et al. multivariate discriminant analyses on plant traits (see below) we tested for multicollinearity among predictor variables using a variance inflation factor (VIF) procedure. Multicollinearity was computed using the vif r package (Lin, 2012) and it allowed us to discard five variables of the initial dataset that showed strong multicollinearity. We set a threshold of VIF = 3, and we discarded redundant variables above this value. Two variables were removed at the plant level (plant diameter and plant cover), and three variables at the cone level (cone mass, pulp mass and total seed mass; Table S2). We tested for differences in cone traits among plants of the three stands by means of linear discriminant analysis carried out with the ‘lda’ function of the mass package (Venables & Ripley, 2002). We used discriminant analysis to find a linear combination of traits that maximizes variation among plants within the three distinct sampled stands. Traits, when combined with others, provided the best discrimination power were selected through the Wilks' lambda criterion, using the ‘greedy. wilks’ function in the R package klar (Heino et al., 2015). 2.6.1 | ERGM theory ERGM has been used recently in ecological research for the analysis of individualbased mutualistic interaction networks (Arroyo- Correa et al., 2021; Miguel et al., 2018). The ERGM design is analogous to a GLM and implements a Markov chain Monte Carlo maximum likelihood parameter estimation. These models allow us to test hypotheses about the underlying mechanisms shaping networks by modelling how the number of links is affected by specific predictor variables associated with the network nodes (node covariates). In common with other models, ERGMs require a link function that captures the underlying distribution of the data (here, interaction strength distribution). In addition to node covariates, to take into account aspects of network topology like link weights, the ERGM allows to consider these structural covariates in parameter estimation. The structural term representing the total number of link weights, 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. 2.6.2 | ERGM interpretation In a weighted ERGM, the coefficient estimates of node covariates must be interpreted in terms of the effect on the probability of increasing (or decreasing) the strength of pairwise interactions. They indicate the expected change in the logged number of any pairwise interaction increasing (or decreasing) as a function of a given oneunit change in a specific node attribute (Arroyo- Correa et al., 2021; Kolaczyk & Csárdi, 2014). In our case these predictor variables are directly associated with the node (plants) attributes (e.g. cone crop size, seed mass) and model how these traits affect the number of links that each node establishes in the network. By associating network topology to specific plant attributes, we aimed to assess the contribution of each node (plant) attribute to the overall network configuration. The underlying Poisson distribution of interactions allows the interpretation of the coefficients as the increase or decrease in the frequency that each plant– animal pair interacts (logged). In this way, plant traits, as node attributes, affect animal preferences that ultimately define the interaction weights measured; these in turn determine how each node of the network (a plant) gets its position in the overall network. 2.6.3 | ERGM design We fitted ERGM with the structure of the individualbased bipartite network as the response variable, the sum of the link weights as a structural covariate, and the plant and cone traits and neighbourhood attributes of individual plants as node covariates. We use the full network to fit the best possible model based on AIC criteria. For this, hierarchical models were built by adding the different groups of variables or factors (area factor, general plant variables and cone quality variables). The best model for the general network is the one that was also used to evaluate the effect of node attributes on the specific networks of each area. Therefore, we finally fitted the same model for the general network and, separately, for the three perstand networks, adding neighbourhood attributes (specific to each area) as an additional node variable. In order to examine the relative importance of each group of variables (neighbourhood attributes, individual plant traits and cone quality traits) in our models, we calculated the combined effect accounted for by each variable group, relative to the summed value of the estimates. The computation of effects variance partitioning and of goodness- of- fit model estimates are not yet implemented for bipartite weighted ERGM (Kolaczyk & Csárdi, 2014), and we did not estimate them. We use the relative sizes of estimates (and the summed estimate for variable groups) as a proxy of the variable relevance in the model (Arroyo- Correa et al., 2021), and to discern between effects due to link configuration and effects given by node attributes. 3 | RESULTS 3.1 | Predispersal seed– pulp predators We found the same assemblage composition of seven predispersal seed– pulp predators species in the three study stands (Figure 2). We quantified their importance through two components of the interaction, frequency of occurrence and prevalence. In these terms, the most prevalent interaction was that of C. chloris with the plants at the colonization front, interacting with 91% of the plants (Figure 3a). However, considering the 105 plants surveyed, the most widespread predators were invertebrate species (Table S3). Overall, in terms of frequency of occurrence of the interaction, C. chloris was the most frequent consumer of J. phoenicea seeds, specifically in the colonization front (Figure 3a, Table S3). We counted a total of 7458 cones in seedtraps, preyed upon by C. chloris, resulting in 1865 estimated interaction 13652745, 2022, 9, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2745.13942 by Universidad De Sevilla, Wiley Online Library on [17/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 2196 | Journal of Ecology ISLA et al. events. Chloris chloris interactions were much more frequent (and also spread between plants) at the colonization front than in the intermediate and mature stands (Figure 3a). From the 7905 collected and dissected cones we found 2175 unique pairwise interactions between juniper cones and invertebrate species. Except for the case of C. chloris at the colonization front, invertebrates were more important in terms of interaction intensity and prevalence (Figure 3a). Finally, rodents (A. sylvaticus and M. spretus) were the least frequent group (Figure 3a; Table S3). We quantified a total of 154 potential rodentfocal plant interaction events, 36% sampled by live traps and 64% by camera traps. 3.2 | Individualbased, antagonistic interaction networks of J. phoenicea The three standbased networks were small and their structure differed from stand to stand (Figure 2). Hamming distances indicated that the networks of mature (MAR) and intermediate (OJI) stands were more similar to each other than compared to the colonization front (COL) network. The minimum number of link substitutions required to change one network into the other was 3760 between COLOJI, 3780 between COLMAR and 2636 for OJIMAR. At the node level, the spread with which the strength of interactions between individual plants was distributed differed between stands. As a result, fewer plants at COL, compared to MAR and OJI, showed antagonistic interaction strength values >0.20 (COL = 8, MAR = 12, OJI = 14, Figure 3b). The node strength of animal species varied greatly between the three networks (Figure 4a). Chloris chloris had a high importance in the COL area, with a reduced importance of the rest of the species in comparison to MAR and OJI, where invertebrate species were dominant. The wasp M. amicorum and the moth M. oxycedrella show a peak in OJI although they were also relevant in MAR. The three estimated network descriptors (wC, wNODF and M) showed marked, almost linear, trends along the landscape gradient (Figure 4b), with significant differences among them. Weighted connectance (wC) increased strongly from MAR to OJI and to COL, when approaching the colonization front, showing a linear trend (Figure 4b, Table S4). The same pattern was also found for weighted nestedness (wNODF) which was significant for the three networks (in comparison with randomly built networks), and also with highly significant differences between the three sites (Figure 4b). In contrast, modularity decreases when approaching the colonization front (Figure 4b). For both wC and wNODF, the MAR and OJI measurements were closer to each other than the COL values. We did not detect any clear pattern FIGURE 2 Weighted bipartite networks of the three stands representing the focal plants of Juniperus phoenicea in each stand (coloured nodes) and the predispersal seed– pulp predator species (black nodes). Each pairwise plant– animal interaction is represented by a grey link whose thickness varies according to its frequency of occurrence. Node labels identify each focal plant. Animal species codes: APO - Apodemus sylvaticus, MUS - Mus spretus, CHL - Chloris chloris, MES - Mesophleps oxycedrella, ELA - Elachista sp., NAN - Nanodiscus transversus, MEG - Megastigmus amicorum. Bottom, illustrative representation of the three stands along the regeneration gradient, from the most mature stand (MAR; a), to the intermediate stand (OJI; b), and ending at the colonization front stand (COL; c). In addition, some differences between stands, such as neighbourhood density, plant size or fecundity, are represented. M029 M036 M058 M065 M087 M136 M144 M176 M207 M208 M215 M219 M304 M315 M381 M457 M460 M518 M520 M531 M541 M581 M648 M664 M694 M702 M749 M820 M826 M843 M847 M849 M956 M986 CHL APO MUS MES ELA NAN MEG O006 O008 O023 O042 O045 O054 O067 O072 O105 O112 O123 O136 O137 O148 O231 O241 O268 O304 O314 O334 O352 O384 O463 O488 O514 O553 O580 O582 O621 O626 O642 O682 O685 O766 O839 CHL APO MUS MES ELA NAN MEG C002 C014 C017 C018 C024 C050 C054 C055 C056 C057 C065 C075 C083 C086 C106 C107 C111 C119 C121 C126 C131 C140 C141 C150 C153 C162 C165 C168 C172 C174 C176 C193 C285 C316 CHL APO MUS MES ELA NAN MEG (a) (b) (c) 13652745, 2022, 9, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2745.13942 by Universidad De Sevilla, Wiley Online Library on [17/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License | 2197 Journal of Ecology ISLA et al. in the number of modules across the colonization gradient (MAR = 3, OJI = 5, COL = 4; Figure S7). The seedpredator M. amicorum and the pulppredator M. oxycedrella belong to the same module in the three networks. In all three networks, C. chloris formed a single module or was accompanied by other vertebrates, but never shared with invertebrate species. All network metrics were significantly different from those calculated by null models; their observed values, and their detailed confidence intervals can be found in (Table S4). 3.3 | Plant characteristics Linear discriminant analysis showed that individual plant and cone traits are good predictors for specific area classification, indicating significant differences in plant characteristics among areas. The two first discriminant functions assign correctly each plant to its stand in 70.47% of the cases. This means that plants within stands tend to be more similar to each other in terms of traits in comparison with plants from the other stands. The significant independent variables according to Wilks' lambda, were plant level variables (plant height and cone crop size) and cone level traits (seed mass and number of seeds per cone; Table S5). The first discriminant function (LD1) accounted for 83.57% of variance, discriminating individual plants of the colonization front (COL) from the rest of individuals (Figure S8). This first function is characterized by the effect of plant height and plant fecundity, with the younger plants from COL being taller and with larger cone crop sizes. Cone traits also discriminated plants from COL, with more seeds per cone and the heaviest seeds. The second discriminant function (LD2) accounted for the remaining 16.43% of trait variance, contributing FIGURE 3 (a) Role of predispersal seed– pulp predator species in terms of frequency of occurrence in our sampling (lines show ±1 SE), and prevalence (percentage of plants with at least one interaction recorded). The isolines (and their values) represent the product of both components and can be used to illustrate the magnitude of the differences between species roles. (b) Rankings of plant individual strengths (plant codes along the x axis) in each study area (most mature stand, MAR, intermediate stand, OJI and colonization front, COL). Dotted lines illustrate a reference threshold set at 0.2 to highlight the number of plants with node strength above or below it in each network. 44 107 184 276 381 503 644 807 1001 1235 1530 1945 0 10 20 30 0306090100 Prevalence (%) Frequency of occurrence (%) MAR 44 107 184 276 381 503 644 807 1001 1235 1530 1945 0306090100 OJI 44 107 184 276 381 503 644 807 1001 1235 1530 1945 0306 09 0100 COL M956 M087 M849 M144 M664 M581 M843 M826 M694 M749 M136 M215 M381 M460 M058 M986 M315 M541 M520 M029 M531 M036 M648 M304 M457 M207 M820 M702 M847 M208 M176 M065 M219 M518 0.0 0.2 0.4 0.6 0.8 1.0 O580 O304 O067 O042 O766 O334 O148 O514 O072 O682 O137 O054 O123 O642 O112 O685 O023 O314 O231 O582 O105 O008 O045 O352 O621 O384 O006 O268 O839 O626 O553 O463 O136 O241 O488 C075 C024 C083 C150 C140 C126 C106 C285 C057 C054 C065 C176 C017 C316 C014 C172 C050 C002 C165 C131 C174 C056 C141 C121 C168 C086 C153 C162 C193 C107 C119 C018 C055 C111 Node strength OJI COL Individual plants Apodemus sylvaticus Chloris chloris Elachista sp. Megastigmus amicorumMesophleps oxycedrellaMus spretus Nanodiscus transversus MAR (a) (b) 13652745, 2022, 9, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2745.13942 by Universidad De Sevilla, Wiley Online Library on [17/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 2198 | Journal of Ecology ISLA et al. slightly to separation between MAR and OJI in terms of the relation between the number of seeds and their perseed weight. According to the second discriminant function, OJI plants had more seeds per cone than the MAR plants with fewer but heavier seeds per cone. Overall, plants of the youngest stand (COL) were the ones that differed the most in their individual traits, being the FIGURE 4 Metrics of the three individualbased interaction networks of individual Juniperus phoenicea and their predispersal seed– pulp predators in the three stands studied; MAR, mature stand; OJI, intermediate maturity stand and COL, the colonization front. (a) Relative importance of animal species in each network, estimated using the species strength metric. (b) Structural metrics at network level in the three stands. 0 5 10 15 20 MAR OJ IC OL Area Area A. sylvaticus M. spretus C. chloris M. oxycedrella Elachista sp. N. transversus M. amicorum 0.250 0.275 0.300 30 35 40 0.24 0.26 0.28 0.30 0.32 0.34 MAROJI COLMAR OJICOLMAR OJI COL Species strengthWt. connectance (wC) Wt. nestedness (wNODF) Modularity (M) (a) (b) General MAR OJI COL Estimate ± SE z value p value Effect Sum 1.525 ± 0.01 82.418 <0.001 General plant traits Plant diameter1−0.063 ± 0.02 −2.717 0.006 − − + Plant height 0.062 ± 0.02 2.580 0.009 0+0 Plant crop size 0.184 ± 0.02 8.165 <0.001 + + + Cone traits Cone diameter 0.109 ± 0.02 3.753 <0.001 0+0 No. seeds/cone 0.165 ± 0.02 6.517 <0.001 +++ Seed viability −0.117 ± 0.02 −5.715 <0.001 −0− Oneseed mass 0.037 ± 0.03 1.330 0.183 −0+ Neighbourhood attributes Juniper neighbourhood crop size − − + Juniper neighbourhood density +−0 TABLE 1 Summary of general Exponential Random Graph Models (ERGMs, all study stands pooled) together with a schematic summary of the results of the ERGM models in each stand (most mature stand, MAR; intermediate stand, OJI; and colonization front, COL). The model estimates the effects of individual features on the weighted network structure. The ‘sum’ effect can be interpreted as the effect of link weight's distribution on the configuration of the observed network. The righthand side of the table shows the direction of the effect: (+) increasing or (−) decreasing the logged number of interactions. If the effect of the variable is significant the symbol appears in bold, and 0 indicates there is no direction of the effect. 13652745, 2022, 9, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2745.13942 by Universidad De Sevilla, Wiley Online Library on [17/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License