Ecology of plant-animal interactions during plant population expansion processes

1 2 CITACIÓN RECOMENDA: Isla-Escudero, J. 2024. Ecología de las interacciones planta-animal durante procesos de expansión poblacional de plantas. Tesis Doctoral. Universidad de Sevilla, España. RECOMMENDED CITATION: Isla-Escudero, J. 2024. Ecology of plant-animal interactions during plant population expansion processes. PhD. Dissertation. Universidad de Sevilla, España. 3 Ecology of plant-animal interactions during plant population expansion processes JORGE ISLA ESCUDERO Tesis Doctoral Universidad de Sevilla, 2024 4 This thesis received financial support by an FPI grant with code PRE2018-085916 from the Spanish Ministry of Science and Innovation, as well as by the projects CGL2017-82847-P and PID2022-136812NB-I00 from Agencia Estatal de Investigación, Ministerio de Ciencia e Innovación, project PAIDI_2020-P20_00736 from Consejería de Economía, Conocimiento, Empresas y Universidad, Junta de Andalucía and by project LIFEWATCH-2019-09-CSIC-13 funded by European Fund for Regional Development (FEDER-EU). 5 Dr. Pedro Jordano y Dr. Miguel Jácome, CERTIFICAN Que los trabajos de investigación incluidos en esta memoria han sido realizados en su totalidad por Jorge Isla Escudero y son aptos para ser presentados ante el tribunal designado para aspirar al grado de Doctora por la Universidad de Sevilla. Y para que así conste, y en cumplimiento de las disposiciones legales vigentes, extendemos el presente certificado a 22 de enero de 2024. Directores: Fdo. Dr. Pedro Jordano Barbudo Fdo. Dr. Miguel Jácome Flores Tutor: Fdo. Dr. Juan Arroyo Marín 6 7 TABLE OF CONTENTS ABSTRACT……………………………………………………………..………………………………………..…09 RESUMEN………………………………………………………………..…..……………………………........…11 CHAPTER 1. General introduction and Objectives……………………………………..………16 CHAPTER 2. Drivers of individual–based antagonistic interaction networks during plant range expansion ……………………………..……………..……………….…..37  Appendix A. Supplementary information from Chapter 2 ……..………..…..69 CHAPTER 3. The turnover of plant–frugivore interactions along plant range expansion: consequences for natural colonization processes …………….......82 – Appendix B. Supplementary information from Chapter 3……..……………..112 CHAPTER 4. Animal-mediated seed dispersal and the demo-genetic configuration across plant colonization gradients …………………..………...…..131 – Appendix C. Supplementary information from Chapter 4……..…………....159 CHAPTER 5. General Discussion, Future Research Lines and Conclusions ...…..174 BIBLIOGRAPHY………………………………..…………………………………………..……..……………185 AGRADECIMIENTOS / ACKNOWLEDGMENTS……………….……………….…….…..…201 8 9 Abstract Plant-animal interactions are crucial for maintaining ecosystem functionality and biodiversity. Ecological research has shown that these interactions can be modified in different contexts, affecting natural dynamics. Currently, significant changes in the landscape, such as rural abandonment or the protection of natural areas, are creating a new ecological context for the formation and regeneration of new forested areas. Although the of the role of mutualistic and antagonistic interactions between plants and animals in ecosystems is acknowledge, we lack a comprehensive understanding of their role during the processes of colonization and expansion. Understanding the structure of ecological interactions and the role of the species involved along natural expansion gradients is key to comprehending the natural dynamics underlying expansion processes, as well as the mechanisms that may or may not facilitate the formation of new forests and biodiverse natural ecosystems. The aim of this PhD thesis is to document the structure and consequences of plant-animal interactions that are related to plant reproductive success along a natural expansion gradient. To achieve this, we analyze the interactions, both antagonistic and mutualistic, of predispersal seed predation and seed dispersal at the individual plant level and the resulting patterns such as seed rain and its maternal composition along the natural expansion gradient experienced by the Juniperus phoenicea subsp. turbinata population in the Doñana National Park since its protection five decades ago. We found that the interaction networks that shape individual plants, both with antagonistic (predispersal seed predators) and mutualistic (seed dispersers) species, were concordantly reconfigured along the population expansion gradient. The structuring at the individual level was influenced by the animals' preferences for certain phenotypic traits of the plants, as well as the neighborhood and maturity context in which they were found. Interactions with most animal species, including mutualists and antagonists, are well established in mature stands; 16 CHAPTER 1 General Introduction and Objectives 17 Global change and land–use change: Plant expansion and diversity Our planet is undergoing a global transformation marked by unprecedented and rapid transformations. The alteration of land use due to human activities is a significant factor directly impacting biodiversity. Human activities such as agriculture, livestock, and infrastructures has resulted in the widespread loss of natural habitats, posing a threat to the persistence of global biodiversity (Newbold et al. 2015). However, both the technification and industrialization of agricultural processes in recent decades have led to an increasing migration of populations from rural areas into urban centers. This phenomenon is particularly evident in more developed regions such as Europe and North America and anticipated to persist in the near future (Johnson & Lichster, 2019, Perpiña Castillo et al. 2018). Rural land abandonment has triggered significant landscape changes. For instance, Europe has witnessed a notable transition, the so-called `forest transition´, shifting from a period of decreasing forest cover each year to a phase of increasing forest cover (Kauppi et al. 2018). Many of this emerging green natural cover did not undergo intentional planting; instead, originated from natural forest regeneration as an expansion of secondary forest (FAO 2020). Understanding the mechanisms determining the success of these expansion processes is crucial for assessing the potential of natural regeneration to colonize new areas and ensure the functionality and diversity of emerging forests. In plant population expansion scenarios, reproductive individuals in remnant areas serve as crucial sources of propagules (seeds) necessary for colonization. Pioneer species, defined as those first colonizing new habitats created by disturbance (Dalling 2008), play a particularly relevant role in this process. These species establish populations that facilitate subsequent successive processes of other species, representing a key step in the formation of mature and stable communities over time. In this context, ecological interactionscrucially those involved in seed dispersal processesare central to shaping the various ways in which expansion processes unfold (Hampe 2011). Such range expansions depend 18 on the populations established at the colonization front, or ’leading edge’ (Thuiller et al. 2008). Hence, assessing and predicting plants expansion potentials requires a sound understanding of the ecological and demographic (demogenetic) processes that occur in such populations (Hampe & Petit 2005, Hampe 2011). These interactions among plants and other elements (e.g., animals) of the ecosystem influence the dynamics of colonization and community dynamics. To gain a comprehensive understanding of these processes, field studies are essential to provide baseline information at local–scales and identify and protect mechanisms supporting natural regeneration. Plant–animal interactions across plant life cycle during plant expansion Plants and animals have coevolved together, leading to a multitude of complex plant/animal interactions (Herrera & Pellmyr 2002, Thompson 2005) involving mutualistic effects. An individual plant establishes multiple interactions with diverse animal species throughout its life cycle (Figure 1.1), being individuals, not species, the biological entities actually interacting in nature. From a phytocentric point of view, plant–animal interactions result in different outcomes on the fitness of individual plants during the sequential stages of each reproductive episode (Harper 1977). Antagonistic interactions, for instance, are those interactions having a negative effect on plant fitness, reducing its reproductive success (Sallabanks & Courtney 1992, Strauss & Irwin 2004). For example, interactions with herbivorous animals that may directly or indirectly affect plant attractiveness to pollinators (Bronstein et al. 2007, Moreirla et al. 2019). Seeds, and the surrounding pulp in fleshy–fruited plants, may also interact with predispersive antagonists, which reduce plant crops before propagules become dispersed (Janzen 1971, Guido & Roques 1996, Crawley 2000, Hulme & Benckman 2002). In addition to antagonisms, most plants benefit from the interactions with animals and the mutualistic relationships they establish with them. For example, interactions with frugivorous animals that feed on fruits facilitate seed germination (Traveset 1998) and sustain seed dispersal to adequate microhabitats (Jordano & Schupp 2000) away from the mother trees, where it is easier to establish 19 and survive (Janzen 1970). These are just a few examples, but there are many other types of plant-animal interactions that are distributed along a gradient from purely mutualistic to fully antagonistic outcomes (Rodríguez–Rodríguez et al. 2017, Jácome–Flores et al. 2020). For every individual plant in nature we can locate the pooled results of the mixing of antagonistic and mutualistic interactions, in terms of fitness effects, along such a gradient (Gómez et al. 2023); that location depending, ultimately, on the way individuals “build” their interactions with animal partners. Figure 1.1. Schematic figure representing some of the main plant–animal interactions that determine the reproductive success of a given plant individual throughout the different stages of its life cycle (see Harper 1977, Wang and Smith 2002). In the field, when we observe and record an interaction (e.g., a butterfly probing a flower), we observe a pair of individuals of different species interacting with each other. In fact, individuals, not species, are the actors of ecological interactions in nature. Using this individual–based approach can avoid underestimating inter–individual variability, which can be more significant within 20 a species than between them (Albert et al. 2010, Soefert et al. 2015). At this hierarchical level, interactions are driven by individual plant context (neighborhood attributes) and individual plant phenotypic traits (Rodríguez- Rodríguez et al. 2017, Bartomeus et al. 2016, Guerra et al. 2017, Snell et al. 2019, Arroyo-Correa et al. 2021). Using the example mentioned earlier, the butterfly may preferentially pollinate individuals growing among heterospecific neighborhoods with high flower synchrony, as opposed to those in monospecific patches (neighborhood attributes); within these patches, it may preferentially visit individuals with higher number of flowers (phenotypic traits). These extrinsic and intrinsic characteristics of individuals vary among populations and actually define the context dependency of most ecological interactions (Chamberlain and Holland 2008). For instance, a plant expansion gradient is a good example of such context variability and can capture this heterogeneity (e.g., neighborhood density, fruit productivity, etc.) (Gratani 2014). Therefore, plant–animal interactions can be expected to be reconfigured along these gradients, either slowing (antagonisms) or accelerating (mutualisms) regeneration, although this topic is still unexplored in ecological research (Svenning et al, 2014). Plant–animal interactions that will have the strongest effect in favoring or limiting the expansion of plant populations are those closely related to their individual propagule contribution (size of the seed crop) and dispersal capability. The loss of propagules which results from pre–dispersal predation by invertebrates, birds and mammals (Hulme & Benckman 2002, Janzen 1971) can trigger a seed source limitation (source limitation, sensu Clark et al. 1999), hindering successive regeneration and expansion stages. Evidently, a local scenario of natural expansion is the result of a successful regeneration process, in which antagonistic agents may not represent a drastic source limitation. One of the possible hypotheses linking natural expansion processes and plant–animal pre–dispersive antagonistic interactions is the enemy release hypothesis (Keane & Crawley, 2002). This release from enemies in stands or populations at the advancing, expansion front (Petit et al. 2004) facilitates the expansion through a reduction or elimination of antagonistic agents, as repeatedly documented for expansions of alien species 21 (Meijer et al. 2016). Following this hypothesis, it is expected that the colonization front stands of expanding plant populations will be exposed to low antagonist pressure, which facilitates local regeneration and expansion beyond the population. However, there is no empirical evidence available to support this hypothesis during plant expansion, nor which are the drivers that shape these critical interactions along population expansion gradients. Moving forward in the plant reproductive cycle, mutualistic seed dispersal interactions are pivotal in plant species expansion (Howe et al. 2004, Hampe 2011, Escribano-Avila et al. 2014). Frugivores provide dispersal services to many plant species in exchange for food resources: 70–94% of woody species in tropical forests and 50–75% in Mediterranean habitats (Jordano 2014). As a result of this interaction, fleshy–fruited plant species are able to colonize new habitats via seed rain generated by frugivores. As with pre–dispersal seed predation interactions, seed dispersal interactions are also driven by both traits of individual plants and the neighborhood context in which they are (Schupp et al. 2019). The combination of interactions with multiple dispersing species that an individual plant establishes will largely determine its fitness, i.e., its contribution in terms of new individuals for population regeneration and expansion to new sites. In addition to the interactions themselves, effects on individual fitness will result from the quantitative and qualitative component of the interaction effectiveness that different frugivores perform (Schupp 1993, Schupp et al. 2017, Quintero et al. 2023). Therefore, to fully understand the role of mutualisms in individual fitness, it is necessary to characterize; i) with whom they interact, ii) how much they interact, and iii) what this interaction results in. The expansion of local plant population ranges is supported by individual contributions to the seed rain (and the success of posterior demographic processes). The release of dispersal–limiting constraints (Nathan & Muller– Landau, 2000, Schupp et al. 2002) can trigger a substantial advance of colonization, especially if the effect of mutualistic interactions exceeds that of antagonists (compensatory effects). However, we do not yet have a good understanding of how 22 such interactions are configured along expansion gradients, nor which individuals contribute more to the seed rain along the gradient. If seed rain comes mainly from individuals of the remnant stands, the expansion process may be somewhat slowed (`wave front expansion´, Skellam 1951). Nevertheless, combinations of short and long dispersal events coming from the whole maturity gradient will accelerate plant population expansion (Bialozit et al. 2006, Excoffier et al. 2009). However, it is important to consider that such accelerated processes may pose a number of important risks to the diversity of future populations. Founder effects coupled with rapid expansion may result in a diversity erosion in the colonization fronts of expanding populations by allowing a few colonizing individuals to establish and reproduce over vacant lands while hampering the arrival of later individuals (Waters et al. 2013). So far, research assessing these questions has mainly been based on mathematical models, with little or no empirical support. In fact, the role of animal–mediated seed dispersal plays in reducing or maintaining the diversity of expanding populations is so far completely unknown. Further understanding of range expansion dynamics can be gained by characterizing the maternal progenies in seed rain, such as the number of source trees reaching specific deposition sites. This approach represents an ideal opportunity in order to assess diversity patterns along expansion gradients (Petit et al. 2004). The maternal composition in the seed rain provides the starting spatial template on which new forests will form (Jordano & Godoy 2002). Relatively few studies have assessed the dispersal of maternal progeny in the seed rain generated by frugivores (Grivet et al. 2005, García et al. 2009, Karubian et al. 2010, Scofield et al. 2012, Carvalho et al. 2020), highlighting the importance of analyzing this aspect by demonstrating the crucial role of frugivorous animals in the maternal composition of the seed rain. For example, the genetic consequences of longdistance seed dispersal events during episodes of range expansion have been explored only recently, mostly using simulations (Nichols and Hewitt, 1994, Ibrahim et al. 1996, Le Corre et al. 1997). The actual link between animal-mediated seed dissemination by frugivores and the maternal composition of the seed rain are the foraging patterns of frugivores (foraging behavior and microhabitat 23 preferences, Figure 1.2). In this way, the demographic and genetic consequences of animal–mediated seed dispersal are tightly intertwined (García and Grivet 2011). Nevertheless, how this maternal composition is structured along plant population expanding gradients, as well as the role of frugivore species (with different foraging patterns) in these ecological contexts, remains an unknown aspect in ecology. Colonization fronts often exhibit sparse fruiting tree populations, potentially diminishing maternal source tree diversity, and thus, the genetic diversity of populations during expansion (founder effects, Waters et al. 2013). However, these effects are not empirically proven and remain hypotheses; for example, even infrequent long-distance dispersal events may suffice to dampen founder effects during early colonization (Bialozyt et al. 2006). If dominated by vectors with generalized foraging patterns (Figure 1.2), mutualistic plant–animal interactions may substantially diminish these negative founder effects, favoring the diversity of the forests of the future. In fact, there are already studies showing that new forests that have formed rapidly may be at least as taxonomically and functionally diverse as mature remnant forests (Espelta et al. 2020, García et al. 2020). The increasing data on forest regeneration and expansion dynamics offer to ecologists a timely opportunity to examine the role of plant–animal interactions during these increasingly relevant processes. 24 Figure 1.2. Schematic figure illustrating contrasting foraging patterns that frugivore species, or individuals, may exhibit. Frugivores transit from fruiting trees, where they forage on fruits, to deposition sites where they rest, perch, or feed again. In this process, they disperse seeds into deposition sites forming strongly aggregated seed rain patterns, typically located at preferred microhabitats (Jordano & Schupp, 2000). Frugivores exhibit different foraging patterns defined by the combination of microhabitat preferences and species-specific foraging behaviors that encompass the sequence of movements among fruiting trees and deposition sites. For example, the blue lines represent a generalized foraging pattern since the frugivore visits a variety of microhabitats (pine trees, fruiting shrubs, open ground) where it feeds and disperse seeds. This foraging pattern maximizes the mix of progeny in space and concentrates it in preferred microhabitats. In contrast, the orange line represents a territorial foraging pattern, in which the bird feeds in a unique shrub, generating a very restricted seed shadow dominated by a single progeny, impoverishing the spatial diversity of the seed rain. 25 A multidisciplinary approach to advance knowledge into the role of plant–animal interactions during expansion processes. The approaches and methodologies that we use to address the study of complex systems largely determine our capacity to understand them and the extent of our scientific inferences. The quantification, characterization and analysis of plant-animal interactions in nature poses a challenge in ecological research homologous to biodiversity measurements. To this end, it is essential that ecologists draw on the knowledge of those who have gone before: “We are like daf  he hlde f gia (Bernardo de Chartres, XII, Chen 2003). However, it is essential to develop new frameworks, methods and tools to monitor, analyze and evaluate the diversified plant-animal interactions in order to move forward without losing sight of the past. An essential characteristic of sampling methods for frugivore-plant interactions is that any method, taken individually, biases to some extent the recording of certain species (e.g., canopy dwellers vs. understory species, large-bodied vs. small species, diurnal vs. nocturnal, etc.). This is so due to the enormous diversity of frugivore species, and a useful approach has been the combination of sampling methods. Sampling plant-animal interactions Sampling methods for plant–animal interactions monitoring are diverse and their application depends on the context, the targeted interaction and the desired data (e.g., interaction outcomes). The development of new techniques in recent decades like molecular ecology tools (e.g. DNA-Barcoding) or multimedia monitoring (phototrapping) has revolutionized the sampling of interactions (González-Varo et al. 2014, Zhu et al. 2021, Arstingstall et al. 2021, Naqvi et al. 2022, Garrido et al. 2023). Furthermore, current advances in Artificial Intelligence and Machine Learning algorithms continue along this path, facilitating many tasks for ecologists (Villalva et al. 2023). These new advances will undoubtedly become key 32 Figure 1.6. Contrasting stands along the plant population maturity gradient. In the top, a detail of the top canopy of the mature stand (Sabinar del Marqués). In the bottom, the stand in the colonization front (Sabinar de Colonización) grows in an area dominated by low monte blanco woody shrubland. Photographs by Pedro Jordano (above) and Jorge Isla (below). 33 Objectives The general objective of the present thesis was motivated by the ecological context, knowledge gaps and approaches presented in the previous section. Our main goal was to evaluate the implications of plant–animal interactions during plant population expansion, i.e., how the complex network of interactions with antagonistic and mutualistic animals is reshaped during colonization. In order to assess this topic, we analyzed two different interactions (pre–dispersal seed predation and seed dispersal) and the demographic and genetic consequences of seed dispersal along a population maturity gradient. The questions addressed here constitute four independent yet interconnected chapters (three research articles and a general discussion) that are arranged sequentially according to the natural progression of ecological interactions in the plant reproductive cycle (Figure 1.1). This thesis represents a step towards our understanding of how plant-animal interactions are shaped in natural plant expansion contexts. However, it is part of a more ambitious goal for ecologists: understanding what patterns govern ecological interactions in nature and what implications these interactions have had, have and will have in a rapidly changing world. Below, I list these chapters, as well as our main objective, hypotheses and the approach used in each of them. Chapter 2. Drivers of individual–based antagonistic interaction networks during plant range expansion. In this chapter we focus on pre–dispersal seed– pulp predation interactions by diverse antagonistic species. We monitored these interactions at the individual–plant level during two years in contrasting stands along a natural expansion gradient. We assess the effect of individual plant traits and neighborhood context as drivers of network topology by ERGM´s network analysis. We expected a reconfiguration of interaction networks along the colonization gradient that would explain the rapid expansion of the population (e.g. through a release of antagonistic interactions in recently established stands, Keane & Crawley, 2002). Furthermore, we hypothesize that the individual–based 34 network topologies are being driven by individual traits and neighborhood context. This study is published in Journal of Ecology (Isla et al. 2022). Chapter 3. The turnover of plant-frugivore interactions along plant range expansion: consequences for natural colonization processes. Here, we focus on mutualistic interactions between individual plants and frugivore species (seed dispersers) along the expansion gradient. We assess the reconfiguration of frugivory interactions along this gradient and the role of intrinsic and extrinsic individual variability in shaping interactions assemblages. Furthermore, we estimated the individual propagule contribution by estimating both the quantitative and qualitative components of plant–animal interactions. We expected to identify a reconfiguration of interactions and/or individual propagule contributions variation along the maturity gradient that can provide clues about the expansion history of a natural population. In addition, we hypothesized that antagonists (Chapter 2) and mutualists (Chapter 3) may have preferences for similar plants (individual traits and neighborhood context) and therefore dispersers may compensate for the negative effects of pre–dispersive antagonists. This study is published in Proceedings of the Royal Society B (Isla et al. 2023). Chapter 4. Animal-mediated seed dispersal and the demo-genetic configuration across plant colonization gradients. Advancing into the regeneration processes, in this chapter we assess the contribution of animal frugivores dispersing seeds across the expansion gradient, from both demographic and genetic perspectives. We identified frugivore species contributing to the seed rain throughout the landscape by DNA-Barcoding and applied genetic markers (SSRs) for maternal analysis of dispersed seeds to investigate how (1) stand maturity, (2) microhabitat types, and (3) frugivore foraging patterns shape the distribution of the maternal progenies along this gradient. As colonization fronts often have sparse fruiting trees, we hypothesize that the negative consequences of founder effects are pronounced in newly established stands. Therefore, we expect that the maternal richness and diversity in the seed rain will be diminished at the colonization front. This study is on the 2nd. review in Journal of Ecology. 35 Chapter 5. General discussion, future research lines and conclusions. This chapter synthesizes and discusses the results of the previous chapters from a broader and more encompassing perspective. It outlines potential future research lines that have been derived from the development of this thesis and details the main conclusions. 36 37 CHAPTER 2 Drivers of individual-based antagonistic interaction networks during plant range expansion Jorge Isla, Miguel Jácome-Flores, Daniel Pareja & Pedro Jordano (2022) Journal of Ecology 110, 2190-2204 DOI: 10.1111/1365-2745.13942 Shortlisted article for Harper Prize Early Career Researcher Award, from the British Ecological Society 38 ABSTRACT 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. Here, we study the changes in individual-based pre-dispersal seed-pulp predator networks along a colonization gradient in a rapidly expanding Juniperus phoenicea population in Doñana National Park (SW Spain). Additionally, we analyzed the role of individual plant traits and neighborhood attributes in network configuration by using Exponential Random Graph Models (ERGM´s). We found seven seed-pulp consumer animal species that varied significantly in their frequencies of interaction and prevalence. While invertebrate species were well established in mature and intermediate stands, the 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. Furthermore, individual plant traits strongly determined the topology of these networks, although with differences between stands. Increasing individual crop size and number of seeds per cone increased the interaction odds of individual plants, while seed viability showed the opposite effect. The topology of the colonization front network 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. Overall, the turnover of antagonistic interactions characterized the colonization front, resulting in more heterogeneous interaction strengths among individual plants. We found no evidence for an antagonistic release of J. phoenicea at the colonization front promoting this rapid expansion. It becomes necessary to explore interactions with seed dispersers to understand how antagonistic and mutualistic 39 plant-animal interactions balance during range expansion. Our study highlights the importance of an individual-based 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, seed-predation. 40 Introduction Fast-paced 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 animal-mediated 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 spatio-temporal 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 41 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). Pre-dispersal 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, pre-dispersal 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 non-legitimate seed dispersers including insects, mammals and birds (Hulme & Benckman 2002, Janzen 1971). The magnitudes of pre-dispersal 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 48 was present in a cone relative to all the analyzed cones, or how frequently we detected at least one C. chloris visit during weekly seed-traps checks. Plant characteristics We thoroughly sampled a set of intrinsic and extrinsic plant traits known to influence animal preferences across the 105 focal plants (Table A.2). We considered plant traits that could drive interactions with pre-dispersal seed-pulp predators hierarchically (Sallabanks 1993): first, the individual neighborhood context, then general plant traits and finally at a higher scale of detail, cone quality traits. To characterize the juniper neighborhood 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 hand-counter. 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, one-seed mass, number of seeds per cone and seed viability was estimated by flotation procedure (Table A.2, Figures A.3-A.5). Individual-based, antagonistic interaction networks of J. phoenicea We constructed four weighted, individual-based, 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 49 individual-based scale is only applied to plants and their individual assemblages. We used the three per-stand 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 plant-level traits in the network topology. In addition, per-stand networks were also used to assess among-stand variation in the role of individual traits and neighborhood attributes on the detected network topologies across the colonization gradient. To detect structural changes in interaction organization among 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 node-level 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 network-level and node-level metrics were computed with R package bipartite, and in the case of 50 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 bipartite package, 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. Assessing drivers of network topology To analyze the effect of individual and neighborhood 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 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 A.2). 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 51 through the Wilks lambda criterion, using the ‘greedy. wilks’ function in the R package klar (Heino et al. 2015). ERGM theory ERGM has been used recently in ecological research for the analysis of individual-based 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. 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 one-unit 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 52 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. ERGM design We fitted ERGM with the structure of the individual-based bipartite network as the response variable, the sum of the link weights as a structural covariate, and the plant and cone traits and neighborhood 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 per-stand networks, adding neighborhood attributes (specific to each area) as an additional node variable. In order to examine the relative importance of each group of variables (neighborhood 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. 53 Results Pre-dispersal seed-pulp predators We found the same assemblage composition of seven pre-dispersal seedpulp predators species in the three study stands (Figure 2.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 2.3a). However, considering the 105 plants surveyed, the most widespread predators were invertebrate species (Table A.3). 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 2.3a, Table A.3). We counted a total of 7458 cones in seed-traps, preyed upon by C. chloris, resulting in 1865 estimated interaction 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 2.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 2.3a). Finally, rodents (A. sylvaticus and M. spretus) were the least frequent group (Figure 2.3a, Table A.3). We quantified a total of 154 potential rodent-focal plant interaction events, 36% sampled by live traps and 64% by camera traps. 54 Figure 2.2. Weighted bipartite networks of the three stands representing the focal plants of Juniperus phoenicea in each stand (coloured nodes) and the pre-dispersal 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 (A), to the intermediate stand (B), and ending at the colonization front stand (C). In addition, some differences between stands, such as neighborhood density, plant size or fecundity, are represented. Individual-based, antagonistic interaction networks of J. phoenicea The three stand-based networks were small and their structure differed from stand to stand (Figure 2.2). Hamming distances indicated that the networks of mature and intermediate stands were more similar to each other than compared to the colonization front network. The minimum number of link substitutions required to change one network into the other was 3760 between colonization and 55 intermediate stands, 3780 between colonization and mature stands and 2636 for intermediate and mature stands. At the node level, the spread with which the strength of interactions between individual plants was distributed differently between stands. As a result, fewer plants at the colonization front stand, compared to mature and intermediate ones, showed antagonistic interaction strength values >0.20 (Colonization = 8, Mature = 12, Intermediate = 14, Figure 2.3b). The node strength of animal species varied greatly between the three networks (Figure 2.4a). Chloris chloris had a high importance in the colonization area, with a reduced importance of the rest of the species in comparison to the mature and intermediate stands, where invertebrate species were dominant. The wasp M. amicorum and the moth M. oxycedrella show a peak in the intermediate stand although they were also relevant in the mature one. 56 Figure 2.3. (A) Role of pre-dispersal 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, intermediate stand, and colonization front). 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. The three estimated network descriptors (wC, wNODF and M) showed marked, almost linear, trends along the landscape gradient (Figure 2.4b), with significant differences among them. Weighted connectance (wC) increased strongly when approaching the colonization front, showing a linear trend (Figure 2.4b, Table A.4). 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 2.4b). In contrast, modularity decreases when approaching the colonization front (Figure 2.4b). For both wC and wNODF, the mature and intermediate stands measurements were closer to each other than the colonization front stand values. We did not detect any clear pattern in the number of modules across the colonization gradient (Mature stand = 3, Intermediate stand = 5, Colonization stand = 4, Figure A.6). The seed-predator 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 A.4) 57 Figure 2.4. Metrics of the three individual-based interaction networks of Juniperus phoenicea individuals and their pre-dispersal seed-pulp predators in the three stands studied. (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. 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 A.5). The first discriminant function (LD1) accounted for 83.57% of variance, discriminating individual plants of the colonization front from the rest of individuals (Figure A.7). This first function is characterized by the effect of plant height and plant fecundity, with the younger plants from the colonization stand being taller and with larger cone crop sizes. Cone traits also discriminated plants 64 plant nodes (Figure A.6), acting as a super-generalized partner for most individual junipers in the colonization front. We also found high values for modularity, an expected result for antagonistic networks (Cagnolo et al. 2011) and also when downscaling to the individual level (Dupont et al. 2014, Tur et al. 2015). The decreasing trend of modularity towards the colonization gradient could be promoted by two processes: (1) the generalization of C. chloris and high prevalence of their interaction reduced the possible modules and their intensity at the colonization front and; (2) plants' traits heterogeneity of mature stands favored the appearance of modules due to a higher specificity of the interactions (Figure A.9). The modularity analysis also reveals a truly checkerspot pattern in the overall antagonistic interaction strength across sites: distinct groups of individual junipers showed up at the three sites in terms of their sharing of interaction partners. These modules never include both vertebrates (especially C. chloris) and invertebrates or, if including just invertebrates, almost invariably they contain a seed predator (M. amicorum) and a pulp consumer (M. oxycedrella). The better representation of invertebrates in the mature and intermediate stand together with the dominance of C. chloris at the colonization front is probably responsible for the greater dissimilarity of the youngest network regarding the two more mature stands, as detected by the higher Hamming distances. Drivers of network topology In the general ERGM, individual crop size was the main driver of network architecture. In seed dispersal, the number of fruits is the best supported and understood individual driver of variation in interactions (Schupp et al. 2019); more seeds are preyed upon as more cones are produced. This production cost could be profitable since the resource is so abundant that a large part of seeds produced can escape from seed predators. This huge cone production strategy could be explained under a predator-satiation hypothesis (Janzen 1971, Linhart et al. 2014), yet framed here at the individual level. We think that the increase in interactions with plant height may be related to the improved visibility and attractiveness of taller trees to 65 some bird species (Perea et al. 2014), and flying invertebrates (moths and wasps) increasing their encounter probabilities during flight and, therefore, increasing predation probabilities upon cones. Ultimately, this might reflect the pervasive effects of size hierarchies in populations of long-lived plants (Weiner & Solbrig 1984) determining also interaction strength hierarchies. Regarding the diameter of the plant, we find a negative effect on the number of node interactions. In our study area these larger diameters are typical of old plants, irregularly branched, with lower fecundity, smaller cones, etc., and therefore less attractive to predators. As the per-stand results indicate, this negative effect is seen in mature and intermediate stands, where this plant physiognomy is prevalent. At the cone scale, both cone diameter and number of seeds increased network interactions, probably because of their greater attractiveness in terms of pulp (pulp eaters) and seed (seed predators) of these individuals. We expected that plants with a higher proportion of filled seeds would establish more interactions because predators actively choose them (Fuentes & Schupp 1998), but we found the opposite pattern, in concordance with previous studies with similar juniper interactions (García 1998). Two hypotheses might explain this result, although as each encompasses different animal groups, both may be valid: (1) for C. chloris, probably they need to peck more cones on more non-viable plants than on plants with more full seeds, where they find the full seeds rapidly and; (2) for the invertebrates, some plant species induce more frequent seed abortion when detecting insect-damaged fruits or indirectly by ovipositing insects introducing fungi that can cause abortion (Sallabanks & Courtney 1992). Our results of interactions with pre-dispersal seed predators drivers are consistent with previous literature (Hulme & Benckman, 2002, Janzen 1971, Kolb et al. 2007, Sallabanks & Courtney, 1992), although their effect on the topology of emergent and complex networks (e.g. combining vertebrates and invertebrates) had not been previously explored. Networks from each stand showed variations in the sign and the significance of predictive traits, in some cases differing between stands. This result is consistent with the findings of Miguel et al. (2018), where the effects of predictors varied between study areas. Here we only detect a consistent trend in network 66 topologies for plant fecundity and number of seeds per cone, both of which had a positive effect on linkage establishment. The combined estimate (summed effects) of the variables studied showed a marked drop in the colonization front, most likely linked to an expected decrease in the neighborhood effect in this newly formed stand. This lower combined estimate for the pooled effects of variables suggests a reduced predictability of interaction build-up in this stand compared to the intermediate and mature stands. We suggest that the more mature remnants had longer ‘interaction histories’ compared to the colonization front, where individual junipers might be farther from stabilizing their interaction modes. At the more mature stands, with probably more than hundreds of years of interaction history, neighborhood effects and plant characteristics drive the likelihood of link formation for individual plants. In contrast, in the colonization front, probably with <40 years of interaction history, only a few plants form a central subset in the network core, with interactions building up according to variability in fecundity and cone traits. These results may evidence the stabilization time of plant-animal interactions in plant expansion scenarios in terms of species, interaction network structures and trait drivers of these interactions. In the present work we did not focus on the consequences of antagonism for individual plant fitness and thus its contribution to limit the expansion processes. The fact that at the colonization front the spread of antagonistic load is found to be concentrated in a reduced subset of plants, together with the lower importance of invertebrates in this stand, may suggest a release of pressure by predispersal antagonists, favoring fast plant range expansion. However, our results evidence a turnover of antagonistic agents, with the widely ranging, granivorous birds being the main driver of topological network trends. Moreover, functionally, this interaction is likely to result in a much higher seed loss at the colonization front. The trend of these juniper populations in recent years, with a rapid expansion (García et al. 2014) does not support the idea that C. chloris might limit expansion either. To fully understand the relationship of plant-animal interactions and rapid expansion it is essential to analyze seed dispersal by frugivores in the system. Our results suggest that compensation by mutualistic seed dispersers at 67 the colonization fronts (e.g. Holthuijzen & Sharik 1985, Escribano-Avila et al. 2014) and/or effective release from antagonists in certain masting years that favor expansion waves (Linhart et al. 2014) may help understanding these rapid processes of plant range expansion. Furthermore, a similar analysis of the drivers of seed dispersal interactions may be key to assess how animal preference processes exert opposing or confluent selection pressures explaining changes in phenotypes and genotypes along these gradients. Future research addressing variations in the structure of seed dispersal networks, their ecological drivers and their functional consequences in range expansion scenarios will be key to a holistic understanding of the role of ecological interactions in plant range shifts. 68 Data, code and review history availability: Zenodo: https://doi. org/10.5281/zenodo.6583228 GitHub: https://github.com/PJordano- Lab/Juniperus_antagonists_2022 Dryad Digital Repository: https://doi.org/10.5061/dryad.0rxwd bs38 (Isla et al. 2022) Publons: https://publons.com/publo n/10.1111/1365- 2745.13942. Acknowledgements Discussions and comments with Blanca Arroyo and Elena Quintero were especially helpful in the development and analysis of the manuscript, along with feedback provided by Ana Benítez and Francisco Rodríguez, and the valuable suggestions from two anonymous reviewers. Iago Ferreiro assisted with the camera-trap data analysis and Isidro Román and Francisco Carro helped us during rodent sampling. Our work was funded by project CGL2017-82847-P, Spanish Ministry of Science and Innovation and grant SUMHAL, LIFEWATCH-2019-09- CSIC-13, POPE 2014-2020, MICINN through the European Regional Development Fund (PJ). J.I.was supported by an FPI grant PRE2018-085916 from the Spanish Ministry of Science and Innovation. D.P. was supported by a JAE-CSIC grant JAE INTRO 2018 from the Spanish Ministry of Science, Innovation and Universities. M.E.J-F. was supported by Investigadores por México-CONACYT (73215). We thank the logistic and facilities support form ICTS-RBD Doñana and the Doñana National Park administrators for onsite access authorizations. Open access funding enabled and organized by ProjektDEAL. 69 Appendix A Drivers of individual-based antagonistic interaction networks during plant range expansion (Chapter 2) TABLE OF CONTENTS Sampling interactions Invertebrates survey. Detailed description of invertebrate interactions survey. Rodents survey. Detailed description of rodent interactions survey. Chloris chloris interactions survey. Detailed description of C. chloris interactions survey. Figures and Tables Figure A.1. Ripe cones of Juniperus phoenicea subsp. turbinata on the plant. Figure A.2. Adult individual of J. phoenicea in the colonization front. A camera trap was installed recording the focal plant. Figure A.3. Measured values of individual neighborhood attributes per plant in each studied stand. Figure A.4. Measured values of individual traits at the general plant level in each stand. Figure A.5. Measured values of individual traits at cone level in each studied stand. Figure A.6. Interaction module composition for the three stands. Figure A.7. Graph showing the position of the focal plants based on the LDA analysis. Figure A.8. Summary of the three ERGM models fitted for the three networks. Figure A.9. Observed values of interaction specificity (H2) for the three stands studied. Table A.1. Taxonomic classification of the pulp/seed consumer species detected. Table A.2. Individual traits measured for each individual plant. Table A.3. Summary of interaction intensities of the antagonistic species. Table A.4. Network level metrics for three stands. Table A.5. Individual plant traits loadings on the first two (LD1, LD2) linear discriminant functions, estimated with the greedy Wilks’ function. 70 Sampling interactions Invertebrates survey At the beginning of the fruiting seasons we quantified invertebrate interactions by harvesting 50 mature cones from branches spread out on each plant canopy. We stored the cones in transparent tuppers, exposed to indirect natural light until analysis. We dissected 7,906 cones, examined individually through a binocular loupe and recorded the interaction by direct evidence (presence of insect pests in pulp and/or seeds) and indirect inferences (imputable marks or remains of species or the presence of species-specific parasitoid and hyperparasitoid species). We identified antagonist species and their marks based on Roques et al. (1984) and from specific signals detected simultaneously with species presence. Rodents survey To assess rodent interactions, we combined two sampling methods: camera-trap survey and field trapping. At the season start we deployed 8 camera traps per stand that were rotated periodically through the focal, sampled plants. We installed the cameras at 1-m height and 2.5-m away, focusing on the ground and the lower sector of the canopy (Figure A.2). Camera traps were active in video mode for 10 days in each focal plant and then we changed to another 8 focal plants in each stand. This ensured that the whole set of plants were sampled every 5 weeks and included at least 80 nights throughout the study. Camera-trap records totalled 99 rodent sightings for a total sampling time of 201600 h for all the three sites combined. For all video records with rodent presence we counted each observation on a focal plant. Although we did not detect mice feeding on seeds or pulp in all the recordings, we considered all the visits to the plant as a potential interaction event. Additionally, we carried out three sampling campaigns of live-trapping during the rodent’s activity peak in order to supplement the interaction records obtained with the camera-based monitoring. We captured rodents for two consecutive nights 71 under the 105 focal plants. Trapping was set with Sherman traps under each focal plant with a bait (dry bread soaked in olive oil). Each trap was checked out at dusk and the captured mice were identified and released. As in camera trap data, we consider the presence of each species in each focal plant as a potential interaction between the rodent species and the individual plant. We combined the field sampling data with video data to impute to each plant the rodent assemblage composition by proportionally assigning the species recorded in traps to the unidentified species in videos. In cases where we had data from video but not from trapping, the most abundant mice species in each plant stand or the species captured in the most immediate plants was imputed. Chloris chloris survey We quantified C. chloris interactions by using seed traps placed under focal plants to collect the consumed cones falling to the ground during visits by greenfinches. Seed traps consisted of plastic trays (40 x 55 cm, 8 cm height) with small holes (1 mm diameter) that allow the drainage of rainwater and covered with wire mesh (1 cm light) to avoid removal by other animals. During the two sampling seasons, we counted weekly (and resetted the traps after cleaning) the number of damaged cones by C. chloris fallen in the seed traps under all the focal plants (35 revisions). We reviewed the camera trap videos containing C. chloris predation events to assess the cone predation-rate per visit. Based on video records of foraging greenfinches during full visits to the plants (mean [±1SD] of 4.5 ± 2.6 cones consumed per visit), we consider as one visit event for C. chloris in an individual plant every four damaged cones found in the seed traps. 72 Figures & Tables Figure A.1. Ripe cones of Juniperus phoenicea subsp. turbinata on the plant. The red-brown tone of the cones indicates that they are ripe, in their second year of development. Of the cones in the foreground, the upper one shows evidence of attack by a microlepidopteran (possibly Mesophleps oxycedrella or Elachista sp.) whose larvae expel part of the cone pulp (similar to sawdust) when feeding. The cone just below, to the left shows a round hole which is possibly the exit of an adult Megastigmus amicorium after development and feeding on the embryo of one of the seeds of the cone. Figure A.2. Adult individual of J. phoenicea in the most open stand of the colonization front. A camera trap was installed recording the focal plant, pointing both to the ground to detect rodents and to the canopy to detect greenfinch visitation. In addition, a tray can be seen underneath the plant, used to quantify the cones that fall after handling by greenfinches. 73 Figure A.3. Measured values of individual neighborhood attributes per plant in each studied stand. In the neighborhood-fecundity plot (right) the groups of dots represent the frequency (count) of plants showing each fecundity score. Neighborhood cone crop sizes were scored as: 2 = 50-1000; 3= 1000-5000; 4= 5000-10000: 5= >10000. Figure A.4. Measured values of individual traits at the general plant level in each studied stand, each point representing one focal plant (n = 105). 80 Table A.5. Individual plant traits loadings on the first two (LD1, LD2) linear discriminant functions, estimated with the greedy Wilks’ function. Significant values (p .) for the Wilks’ lambda appear in bold. 81 82 CHAPTER 3 The turnover of plant–frugivore interactions along plant range expansion: consequences for natural colonization processes Jorge Isla, Miguel Jácome–Flores, Juan M. Arroyo & Pedro Jordano (2023) Proceedings of the Royal Society B 290, 20222547 DOI: 10.1098/rspb.2022.2547 83 ABSTRACT Plant–animal mutualisms such as seed dispersal are key interactions for sustaining plant range shifts. It remains elusive whether the organization of interactions with seed dispersers is reconfigured along the expansion landscape template and, if so, whether its effects accelerate or slow colonization. Here we analyze plant– frugivore interactions in a scenario of rapid population expansion of a Mediterranean juniper. We combined network analyses with field surveys, sampling interactions between individual plants and frugivores by DNA–barcoding and phototrapping over two seasons. We assess the role of intrinsic and extrinsic intraspecific variability in shaping interactions and we estimate the individual plant contributions to the seed rain. The whole interaction network was highly structured, with a distinct set of modules including individual plants and frugivore species arranged concordantly along the expansion gradient. The modular configuration was partially shaped by individual neighborhood context (density and fecundity) and phenotypic traits (cone size). Interaction reconfiguration resulted in a higher and more uneven propagule contribution, with most effective dispersers having a prominent role at the colonization front stand, where a distinct subset of early arriving plants dominated the seed rain. Our study offers new insights into the key role of mutualistic interactions in colonization scenarios by promoting fast plant expansion processes. KEYWORDS: Mutualistic interactions, plant range shifts, seed dispersal, ecological networks, reproductive ecology, frugivory. 84 Introduction Plant–animal interactions are the cornerstone of many ecosystem functions, crucially determining the persistence of biodiversity (Bascompte & Jordano 2017, Kissling & Schleuning 2015). These interactions are relentlessly shaped and reconfigured by the ecological context in which they occur, leading to variations in their outcomes (Bronstein 2015). For example, interactions may vary in terms of likelihood of occurrence, frequency, sign and strength due to factors such as climate change (O’Gorman et al. , Fricke et al. 2022), elevational range (Plowman et al. 2017), landscape fragmentation (Larsen et al. 2005) and during plant range shifts (Isla et al. 2022). Of strong interest is the role of interactions in plant range shifts, given the rapid changes that natural systems are currently facing (García et al. 2017). In such scenarios the ability of plant populations to provide propagules is essential to withstand and colonize new habitats (González–Varo et al. 2021). Approximately 70–90% of known woody species depend on vertebrate animals for seed dispersal (Jordano 2014), highlighting the critical role of plant– frugivore interactions in colonization processes and plant range shifts (Urban & Keitt 2001). Range shifts occur due to the combined and consistent responses of local populations to a given global change driver. However, the reconfiguration of mutualistic interactions at local scales during the expansion process remains unclear. Recent studies demonstrate that antagonistic plant–animal interactions can be reshaped during local natural expansion scenarios, and that focusing on this scale is important to bridge interaction shifts with their demographic effects (Isla et al. 2022). Studying the reconfiguration of interactions with seed dispersers at the population level over small scales may be key to understanding how this process occurs, providing an empirical link between shifts in interaction network topologies, with the demographic and vegetation dynamics occurring during range shifts. 85 Ecological networks have proven to be an effective theoretical framework and analytical tool for the study of biodiversity architecture (Bascompte et al. 2006). While ecological networks have traditionally focused on species–species interactions, other organizational hierarchies may be more appropriate for addressing some ecological questions (Guimarães 2020). In fact, individuals, not species, are the actors in natural interactions, and therefore, a sub–specific approach can avoid underestimating inter–individual variability, which can be more significant within a species than between species (Albert et al. 2010). Additionally, the individual scale allows for an effective link between network structure and demographic effects. In seed dispersal studies with a phytocentric perspective, individual–based bipartite networks can consist of two different groups (modes) of nodes. A particular plant species is typically represented by individual plants, whereas animals can be represented either by individuals (individual–individual interactions) or species (individual–species interactions). Previous research has shown that this individual scale is a powerful approach for examining various issues in seed dispersal interactions, such as intraspecific specialization–generalization processes (Vissoto et al. 2022, Guerra et al. 2017, Crestani et al. 2019), fitness consequences (Guerra et al. 2017, Jácome–Flores et al. 2020) and network configurations (Friedemann et al. 2022). Intraspecific variability in interactions between plants and seed dispersers can arise from various intrinsic and extrinsic drivers. Plant phenotypic traits are the source of intrinsic variation, with individual fecundity (crop size) and fruit traits being the main drivers of variation (Schupp et al. 2019). Regarding the extrinsic context, variation in neighborhood composition is also a strong driver of intraspecific differences in seed dispersal (Schupp et al. 2019). An ecological context, such as a natural shift in plant range, represents a source of inter–individual variability resulting in variation in plant stand characteristics, individual density, spatial coverage or demographic and genetic structure along the gradient between mature stands and the colonization front (Excoffier et al. 2009, Petit 2011). This intraspecific variability can drive topological patterns of interactions that have the potential to reshape the interaction 86 networks. For instance, modularity is a pervasive property of interaction networks that could reflect habitat heterogeneity (Olesen et al. 2007). In an individual–based seed dispersal network, modules correspond to groups of individual plants sharing a higher number of interactions with similar frugivore species when compared to other plant groups. Recent analyses of modularity in individual–based seed dispersal networks revealed clusters of closely linked plant–animal nodes consistent with the ecological context in which the plants and animals are found. For example, Vissoto et al. (2022) documented distinct subsets of individual plants forming modules with frugivorous species, based on their use of different habitats by the animals. In a similar approach, Friedemann et al. (2022) found that, along a gradient of rainforest habitats, the modular clustering between frugivores and individual palms was largely due to habitat context rather than individual plant traits. Whether the organization of mutualistic interactions during natural plant expansion follows consistent, generalizable patterns remains unknown (Svenning et al. 2014). Assessing mutualistic interaction turnover and its functional consequences, such as propagule contribution, along plant range shifts will help us understand how plant expansion processes occur in nature. In the context of plant–seed disperser interactions, the main outcome from the plant’s perspective is seed movement from the mother plant to other locations. To accurately estimate the number of seeds actually dispersed from an individual plant, a comprehensive understanding of the factors that affect this process is crucial. Details of the interaction such as characteristics of the species–species visit and the seed treatment during ingestion or digestion are essential to obtaining reliable estimates of seed dispersal (Jordano & Schupp 2000). Additionally, an individual–based approach that incorporates data on the interaction between individual plants and frugivores, as well as plant phenotypic traits such as seed viability or the number of seeds per fruit, is necessary for estimating the individuals’ seed contributions (Vissoto et al. ). In natural settings, individual variation in seed contribution to the pooled seed rain is primarily driven by differences in the composition of the frugivore assemblages among plants (Herrera & Jordano 1981). For example, plants that attract diverse frugivore assemblages may 87 have an advantage (Schupp 1993), although a small and effective assemblage may also suffice for successful dispersal (García 2016). Reconfigurations of individual assemblages can have significant implications for seed dispersal and the contribution that individual plants can make to the process of plant–range expansion. In turn, the variability in these contributions can provide insight into how fruit resource selection by frugivores drives the processes underlying range expansion in plants. For example, in situations where individual contributions are very similar, as determined by high similarity in frugivore assemblages, a homogeneous expansion may occur, leading to a continuous spread by diffusion and resulting in a ‘wave front’ (Skellam ). Alternatively, when the action of frugivores is uneven, only a small and distinct subset of individual plants may contribute most of the seed rain. Such strongly unequal contributions of propagules during pollination and seed dispersal have been repeatedly documented for plant populations (Jácome–Flores et al. 2020, Bacles & Ennos 2008, Gerzabek et al. 2020). During a local expansion process, this unequal distribution of individual contributions may be evidence of long–distance dispersal events (Le Corre et al. 1997). These strongly dispersing individuals may serve as pioneer plants that colonized new areas, and contribute disproportionately to local dispersal, and colonization (Bialozyt et al. 2006, Hamrick & Trapnell 2011). By tracking individual seed contribution and disperser assemblages along expansion gradients, it is possible to infer how these expansion processes are occurring. In this study we assess how an individual–based plant–animal interaction network is re–structured by intrinsic and extrinsic factors along a plant expansion process. We also explore the functional consequences of recorded interactions in terms of individual seed contributions. In addition to the above objectives, we believe that analyzing patterns of interaction in these contexts can provide clues about the expansion history of natural populations. Our study system consists of a Juniperus phoenicea subsp. turbinata population and its assemblage of seed– dispersing birds and mammals in the Doñana National Park (Spain). Species of the genus Juniperus are considered as foundation species in Mediterranean ecosystems. The protection of this area five decades ago has facilitated the rapid 88 expansion of this species, generating a natural gradient of plant range expansion from mature stands to colonization fronts in the surrounding shrublands (García et al. 2014). Recent research in the same system indicates that this rapid expansion has not been facilitated by an antagonistic release (Isla et al. 2022), suggesting that seed dispersers may have reconfigured their interactions across these landscapes, promoting rapid expansion from the colonization fronts. This ecological context, combined with the diverse assemblage of birds and mammals known for similar juniper species (Jordano 1993, Escribano–Ávila et al. 2014, Santos et al. 1999), provides an optimal study system to address the following objectives: (i) to evaluate the individual–based frugivory network of J. phoenicea configuration along a plant range expansion gradient; (ii) to explore the role of intraspecific variability and neighborhood context in the observed network structure; (iii) to estimate individual contributions to seed dispersal along the gradient and its drivers; and finally (iv) to explore how this current picture of individual–based interactions along the gradient may provide us with clues about the expansion history. Materials and methods Study site and study species Juniperus phoenicea subsp. turbinata (Cupressaceae) is a large shrub inhabiting coastal dunes and rocky habitats in the western Mediterranean (Adams 2011). Adult trees have an extended fruiting period (September–May) and produce fleshy brown–red cones, also called galbules (Figure 1.4). Juniperus phoenicea cones average five seeds per cone (s.d. = 1.17, range 1–10) and low seed–viability rates (mean ±1 s.d.) of 32 ± 14% viable seeds. Cones act functionally as fleshy fruits and are a key resource for several wintering migrant species like song thrush (Turdus philomelos) (Jordano 1993) and generalist mammals such as red fox (Vulpes vulpes, Figure B.1) and badgers (Meles meles) (Perea et al. 2012). 89 The study was carried out in the Doñana National Park in SW Spain, Figure 1.5). We selected three 1 ha stands along the juniper natural regeneration gradient. All the individual junipers were identified and georeferenced. The mature stand, named Sabia del Ma, has the highest juniper density and is an old forest with a high dominance of J. phoenicea (882 individuals/ha). Advancing towards the colonization front, the second stand called Sabia del Ojill represents a stage of intermediate maturity with a lower juniper density (700 individuals/ha) and composed of Mediterranean scrubland and scattered pine trees. The last stand is named Sabia de Cliaci, and represents a colonization front or expansion area, which is characterized by a low density of scattered individuals (126 individuals/ha). Besides variation in juniper density and intraspecific composition among stands, there are also deterministic differences in individual phenotypes along the maturity gradient. Plants increase in size (height and cover), fecundity (crop size) and cone size toward the colonization front (Figure B.2). In order to select individual plants uniformly distributed within each stand, we established five evenly distributed subplots in each stand (mean ± s.d.) of 522 ± 115 m2 subplot area. We used these subplots to randomly select 35 focal individual plants in each stand (N= 105). This protocol allowed us a stratified–random selection of the focal plants, ensuring a homogeneous distribution of individuals throughout the three stands. Fieldwork was carried out during two consecutive fruiting seasons (October 2018–May 2019 and September 2019–May 2020). Frugivore interactions survey Interactions between focal plants and seed dispersers were sampled using two complementary methods. We combined camera–trap monitoring of animal visits with DNA–barcoding analysis of fecal samples collected beneath the plants. The combination of sampling techniques is especially useful for monitoring diverse assemblages where each method is more adequate for specific interactions (Quintero et al. 2021). 96 Table 3.1. Frugivore species recorded visiting and feeding on cones of Juniperus phoenicea. The table shows the network module assignment, the visits detected in each of the stands, the total number of visits and the consumption rate per visit. The estimated number of seeds dispersed per frugivore and their relative contribution to the seed rain are also shown. Species within each module are ranked by the estimated number of dispersed seeds (decreasingly). 97 Module composition The individual–based frugivory network structure (an individual–species weighted bipartite network) was significantly modular (p–value < 1 × 105), with three densely connected groups of individual plants and animal species (Figure 3.1, Table 3.1). The stand membership of the plants within the modules differed from a random composition (Pearson’s X2 = 22.12, d.f. =4, p–value = 0.0001), revealing a modular organization, highly concordant with the expansion gradient. That is, some modules were significantly over– or under– represented by plants from certain stands (e.g. plants in the colonization front were over–represented in module A, but module C was under–represented). There was a dominant module including most plants from the colonization front (A), a heterogeneous one with a mixed composition (B), and a third module associated with plants from the two most mature stands (C, Figure 3.1). The first module (A) consisted of 39 focal plants linked to T. philomelos and G. genetta. Turdus philomelos was the main frugivore, which accounts for virtually all interactions in the module. This module was found to be significantly overrepresented by plants from the colonization front (Colonization stand = 22, Intermediate stand = 10, Mature stand = 7). The second module (B) was the largest module, formed by 42 individual junipers and three secondary frugivore species in the network (V. vulpes, T. merula and E. rubecula). Module B included a similar number of plants from the three stands (Colonization stand = 13, Intermediate stand = 13, Mature stand = 16). The last module (C) was the smallest one and included 24 plants and the less frequent frugivorous species (T. iliacus, T. torquatus, C. melanocephala, S. atricapilla, O. cuniculus, C. cooki and M. meles). In this module, plants from the colonization front were found to be significantly under–represented (Colonization stand = 0, Intermediate stand = 12, Mature stand = 12). 98 Figure 3.1. Individual–based frugivory network of Juniperus phoenicea along a range expansion gradient. Dots represent individual plants and squares represent frugivore species. As a spatially explicit network, the location of the plants within the panel for each stand corresponds to their actual georeferenced positions. The dotted lines represent a fictitious separation, as the stands are distant from each other (Figure 1.5). Frugivore visits to the focal plants (network links) are represented by grey lines. The colors of the network nodes represent the modules to which the individual plants and animals belong (A = red; B = salmon; C = grey). The three rectangles illustrate the three different study stands along the expansion gradient, from a mature forest area to the colonization front. 99 Plant traits in modules We found that individual phenotype and especially neighborhood attributes could partially predict the module membership of individual plants (Figure B.4). In the LDA analyses, the two linear discriminant functions with the highest discriminating power correctly assigned each plant to the module of origin in 62.85% of the cases. These discriminant functions included two neighborhood– related factors (neighborhood density and fecundity) and one cone trait (cone diameter, Table B.4). By computing these analyses in parallel with the subset of 82 plants with available pulp nutritional data, we found that the most variable chemical trait (pulp’s ash content (), Figure B.) was also a useful trait for correctly assigning the module provenance (Table B.5). Although some traits were effective in assigning plants to each module, the degree of overlap between module and phenotypes was important (Figure B.4). For example, high neighborhood density together with high values of ash concentration differentiated the module C plants from those of modules A and B. Neighbourhood crop size and cone diameter acted in the opposite direction, separating some of the plants from module A (with highly productive neighborhoods and bigger cones) from plants of modules B and C. Finally, the plants in module C occupied a much smaller phenotypic space in the LDA plot, which indicates that the plants from intermediate and mature stands (which compose this module) that interact with this subset of animals share similar characteristics. Functional consequences: individual seed contribution We estimated a total of 884697 dispersed (viable) seeds from the 105 focal plants during the two sampling seasons. Turdus philomelos, which consistently showed a marked role as a connector (within and among modules) of the interaction network (Figure B.6), accounted for 82% of the dispersed seeds. The remaining 18% of seeds were mainly dispersed by T. merula (6,3%), O. cuniculus (5,3%), E. rubecula (3,5%) and V. vulpes (2,3%). The number of seeds estimated to be successfully dispersed by individual junipers differed significantly between 100 modules (linear model; R2 = 0.07, F = 4.72, d.f. = (2, 102), p = 0.01, Figure 3.2). Plants in module A contributed more to the seed rain than plants from modules B and C (Mean ± SE, A= 13 441 ± 2098; B = 5888 ± 859; and C= 4718 ± 628 seeds). Comparing between pairs of modules, the Tukey test revealed significant differences between seed rain contributions of plants from module A and B (pA–B = 0.01, pA–C= 0.08, pB–C= 0.93). We also detected significant differences in the estimated number of dispersed seeds between stands, increasing from the more mature stands to the youngest colonization front area (mean ± s.e., Mature stand = 2724 ± 372; Intermediate stand = 7285 ± 1066; Colonization stand = 15 268 ± 2361; Linear Model; R2 = 0.24, F= 7.22, d.f.  (, ), p .). Comparing between pairs of stands, a Tukey test revealed that the seed contribution by the mature stand plants was significantly lower than those of intermediate and colonization stands (both with p < 0.0001), with no significant differences between intermediate and colonization stands plants (p = 0.42). Figure 3.2. Estimated number of seeds dispersed by each individual plant in the three network modules. Each dot represents an individual plant and the colors illustrate the stand provenance. Note the log scale on the y–axes to improve visualization. 101 Frugivore assemblage composition strongly determined the number of dispersed seeds (Table 3.2). The fitted model accounted for 75% of the variation in individual seed contribution. Turdus philomelos visits had a positive and significant effect, explaining 29% of the variance in the model. Plant reproductive traits also had significant effects, especially individual seed viability, which accounted for more than 17% of the variance in the model. Although less important in the model, visits by small–sized birds and mammals also had a positive and significant effect on individual plant contribution to the seed rain. We found a markedly uneven contribution of juniper individuals to the pooled estimate of dispersed seeds within stands (Figure 3.3). Consistently among stands, few individuals deviated markedly, by several standard deviations, above the stand mean contribution of dispersed seeds Figure 3.3). However, in the colonization front, this ranked distribution of individual contributions was markedly truncated, with very few individuals having a disproportionately high contribution to the pooled seed rain: three individuals in the colonization stand accounted for almost 50% of the seeds dispersed. By contrast, the individual contributions in the intermediate and mature stands were more even, with nine and six individual junipers contributing greater than 50%, respectively. Concordantly, these three individuals occupied a marked central role in the network topology (Figure B.6), with an important function both within their module and also between modules (Figure B.7). We explored the relationship between the contribution to pooled seed rain and the size of the individuals. Using the individual cover as a proxy of plant– age, we detected a positive and significant relationship between individual cover and its seed contribution (linear model; R2 = 0.60, F = 155.80, d.f. = (1, 103), p =2.2  ; Figure .b), revealing a strong size effect. 102 Table 3.2. Multiple linear regression model summary showing the effects of frugivore functional groups visits and plant reproductive traits on individual seed contributions. In addition to the linear regression coefficients for each variable, t–values, p–values and the percentage of explained variance are also shown. 103 Figure 3.3. Ranked z–scores representing the deviation of each individual juniper from the stand mean number of seeds dispersed, i.e. how each individual contribution to the total seed rain deviates from the stand mean along the natural expansion gradient. Points represent individual plants and colors indicate network module assignment. Note that a more homogeneous normal distribution in the contribution would be represented with half of the individuals above the mean and half below the mean (Intermediate stand case), while a strongly skewed distribution with some outliers would result in a strongly truncated line (Colonization stand case). 104 Figure 3.4. A) Comparison between 1977 and 2020 aerial images of the studied colonization front stand. The three individual plants marked are the individuals with an outstanding contribution to the seed rain. This contrast of images demonstrates that these now highly fecund individuals were pioneer individuals already established at this colonization front in 1977. We have validated the correct individual assignment by overlaying geographical coordinates, together with in situ spatial triangulation based on fixed landmarks. The historical image of 1977 comes from plane flights called `Vuelo Interministerial´ carried out by the Instituto Geográfico Nacional (available at: https://fototeca.cnig.es/fototeca/). The 2020 satellite image is available in Google Earth Pro (see https://www.google.com/intl/es/earth/versions/#earth–pro). B) Relationship between the plant canopy cover and the number of estimated dispersed seeds for each plant. The dark grey line represents the fitted log–log linear regression and the grey area spans the slope confidence interval at 95%. Note the log scale on the y–axes to improve visualization. Individual canopy cover here is used as a proxy for plant age. Plant individuals identified in A are indicated. 105 Discussion Frugivore–mediated seed dispersal plays a crucial role in the expansion of plant ranges (García et al. 2017), resulting in demographic consequences at both local and regional scales (Nathan & Muller–Landau 2000). The varying outcomes of frugivore interactions with different species can be linked to demographic effects, which in turn can lead to population–level consequences, including range expansion or retraction. Our empirical findings demonstrate a turnover in frugivory interactions across a local gradient of natural expansion, with consequential effects on variable individual contributions to the seed rain. The entire interaction network was highly structured, featuring a distinct set of modules composed of individual plants that matched their stand provenance and were closely associated with the natural expansion gradient. This indicates that during range expansion events, plant–frugivore interactions are not only shaped by stand characteristics, but they also influence distinct outcomes of the seed dispersal process, such as variable levels of individual plant contributions to range expansion. Our results underscore the importance and interdependence of this mutualism in promoting natural population expansion. Juniperus phoenicea frugivore assemblage The seed disperser assemblage of J. phoenicea was dominated by thrush species and complemented with small birds and mammals. Our results are consistent with those reported previously for other Mediterranean junipers (Jordano 1993, Escribano–Ávila et al. 2014, Santos et al. 1999). However, the extensive sampling effort of this study and the combination of molecular and camera–trap techniques allowed us to describe this assemblage at the individual plant scale. For instance, we were able to detect that Erithacus rubecula had a significant importance (3.5% of dispersed seeds and 11.5% of total recorded visits), despite the weak trait–matching between cone size (mean [range], cone diameter = 9.3 mm [6.2–12.6]) and the E. rubecula gape size (7.8 mm). It may be that E. 112 Appendix B The turnover of plant-frugivore interactions along plant range expansion: consequences for natural colonization processes (Chapter 3) TABLE OF CONTENTS Dna-Barcoding. Detailed description of laboratory procedures for frugivore identification by DNA-Barcoding. Parallel LDA analyses. Detailed description parallel LDA analyses Figures and Tables Figure B.1. Visits of T. torquatus and V. vulpes to J. phoenicea. Figure B.2. Variability of individual neighborhood context and phenotypic traits. Figure B.3. Schematic example of individual seed contribution estimation. Figure B.4. Plants position according to the LDA analysis. Figure B.5. Variability of individual pulp chemical compounds for the 105 focal junipers. Figure B.6. Relationship between seed contribution per individual and the individual centrality. Figure B.7. Network topological roles of individual plants and the 12 frugivore species. Figure B.8. Individual seed viability in each studied stand. Figure B.9. Aerial images illustrating the expansion process of Juniperus phoenicea. Figure B.10. Plants position according to the LDA analysis with chemical nutritional data. Table B.1. Sampling effort specifications at the individual plant level. Table B.2. Individual traits measured for each individual plant. Table B.3. Correction factors used for individual seed contribution estimation. Table B.4. Neighborhood context and individual plant traits loadings LDA with 105 plants. Table B.5. Neighborhood context and individual plant traits loadings LDA´s with 82 plants. 113 DNA-Barcoding We used DNA-barcoding to identify the bird and mammal species that visit focal plants by collecting scats (or regurgitated seeds) in seed traps under the 105 individual junipers. Seed traps consisted in plastic trays covered with 1 cm mesh wire to avoid seed predation (González-Varó et al. 2014). One tray per plant was placed, except for large plants, where two trays were installed to increase the area sampled. Additionally, we delimited a rectangular soil surface next to the tray where we also collected scats to increase the sampled surface. Scats collection was carried out every ten days. Samples were stored at -20ºC and processed following protocols described in González-Varo et al. (2014). Animal DNA was extracted from the surface of the scat or regurgitated seeds, allowing the identification of the frugivore species. In the case of mammal samples, visual identification in the field was also cross-checked with DNA-barcoding identification. Frugivore species identification was based on a 272-bp mitochondrial DNA region (COI: cytochrome c oxidase subunit I). All bird samples were amplified by PCR using the COI-fsd- degF and COI-fsdR primers (González-Varó et al. 2017). Mammals samples were amplified following protocols and primers of Alcaide et al. (2009). PCR product was sequenced and verified for its matching with COI sequences from Barcode of Life Data (BOLD) and Nucleotide Basic Local Alignment Search Tool (BLAST from NCBI) databases. For those samples without successful amplification we performed nested PCR where we used the same primers on the amplicon of AWCintF2/AWCintR4 (Avian DNA barcodes, (Lijtmaer et al. 2012). We considered as a successful identification, all the sequences that showed a >98% percentage of similarity and at least 100-bp. 114 Parallel LDA analyses The subset of data for which the full set of traits was available consisted of 82 focal plants (Colonization stand = 19, Intermediate stand = 32, Mature stand = 31). The predictor variables used for these analyses included two at the neighborhood scale (density and neighborhood productivity), four at the individual plant scale (crop-size, height, diameter and cover), two cone physical traits (diameter and number of seeds/cone), seed viability, and finally four pulp nutritional traits (protein, fibre, total fat and ash content). Since linear discriminant analysis (LDA) does not allow for missing values, we performed a parallel analysis to assess the role of pulp nutritional content. We computed two LDAs with this 82-plant dataset, one including and one excluding nutritional information to assess whether the discriminatory power increased with the addition of these traits. By contrasting the discriminatory power of both LDA's, we found that the analysis with nutritional information correctly predicted module assignment in a higher percentage of plants (71% vs. 65%, Figure B.10). In both models, neighborhood productivity was significant in discriminating between modules. In addition, in the model that included nutritional information, ash content was significant for module discrimination (Table B.4). The parallel LDA analyses allowed us to detect an effect of the ash content. If we would have the nutritional information for the complete LDA (105 plants), we would possibly have found this effect of ash content significant as well. The load of ash content in the discriminant functions (Table B.5) indicates that this trait is especially useful for discrimination between plants from module C and plants from modules A and B. 115 Figure B.1. Top, a group of ring ouzels (Turdus torquatus) perched on a Juniperus phoenicea plant. Bottom, a camera-trap video capture of a red fox (Vulpes vulpes) feeding cones in a focal plant. 116 Figure B.2. Variability of individual neighborhood context and phenotypic traits (at the plant and cone scales) for the 105 focal junipers (coloured dots) in the three stands. 117 Figure B.3. Factors included in the estimation of the individual seed contribution to the pooled seed rain. In addition to the main data, which are the number of visits (recorded by phototrapping and DNA-barcoding), we included information on frugivore foraging (number of cones per visit and percentage of undamaged seeds) and characteristics of each individual plant (the mean number of seeds per cone and the viability of their seeds). This detailed information allowed us to estimate the number of viable seeds dispersed by each frugivore species on each individual plant (right shaded-area), and therefore how much each plant contributed to the pooled seed rain (SEED DISPERSAL CONTRIBUTIONC106). As an example, the thickness of the grey lines represents the relative change in relevance of different frugivore species at plant C106 when considering just the number of speciesplant visits (links to the left-hand side), and the resulting estimate of seeds dispersed by each frugivore species on that individual plant (right-hand side) after accounting for feeding rates, seed viability, and seed damage while handling and digesting cones. 118 Figure B.4. Plants position according to the two discriminant functions computed for the 105 focal Juniperus phoenicea plants, based on plant characteristics and neighborhood descriptors. Each point represents an individual focal plant and colors illustrate the network module to which each plant was assigned. Figure B.5. Variability of individual pulp chemical compounds (% relative to dry mass of pulp) for the 105 focal junipers (coloured dots) in the three stands. 119 Figure B.7. Network topological roles of individual plants (green dots) and the 12 recorded frugivore species (red dots). Following Olesen et al. (2007), the connectivity among modules (c) measures how connected an individual is to other modules and withinmodule degree (z) represents how well-connected an individual is to the individuals within its module. The threshold lines represent the 90th percentile of the observed values in the network for c and z values. The three individuals of the colonization front contributing outstandingly to the seed rain are identified within the dotted ellipse and occupy positions above the population mean in terms of their connectivity within and between modules. 120 Figure B.6. Log-log relationship between the estimated seed contribution per individual to the pooled seed rain and the individual centrality (weighted closeness) in the interaction network (note log scale in both axes). The positive relationship indicates that the most central individuals in the network (with a closer connection to the other frugivore nodes than other plants) are the individuals that make the largest contributions to seed rain. Note that the three plants identified (dotted ellipse), are the individuals accounting for the largest contribution to the seed rain in the colonization front and identified as being already established by 1970 (Figure 3.4) Figure B.8. Measured values of individual seed viability in each studied stand. Note that the three super-fecund individuals of the colonization front are identified to emphasize that the outstanding contribution of these individuals is in part related to their high seed viability. Boxes span the interquartile (25-75%) range, horizontal lines indicate the median value, and lines span the full range of values across individual plants (dots 121 Figure B.9. Aerial images illustrating the expansion process of Juniperus phoenicea at the colonisation front stand in Doñana Biological Reserve. The 1956-1957 and 1977 images are from aerial photographs taken by plane flights, while the 2005 and 2020 images are satellite images. The older images have been obtained from Instituto Geográfico Nacional (https://fototeca.cnig.es/fototeca/). The most recent images from Google Earth pro. Coordinates (36°59'59.80"N, 6°30'23.16"O). The inset rectangle represents the colonization stand, where the 35 focal junipers are currently established. The resolution of the 1950s photo is not sufficient to appreciate whether juniper trees were already present in the area. However, several pioneer individuals are identified in the 1977 image (see Figure 3.4a). The rapid expansion process becomes very clear when comparing juniper density and cover in the area between different years. 128 Table B.4. Neighborhood context and individual plant traits loadings on the first two linear discriminant functions (LD, LD), estimated with the greedy Wilks’ function. Significant values (p  .) for the Wilks’ lambda appear in bold and indicate those variables that are significantly valuable for classifying plants according to the network module which they belong. 129 Table B.5. Neighborhood context and individual plant traits variable loadings on the first two (LD, LD) linear discriminant functions, estimated with the greedy Wilks’ function for the 82 plants including pulp nutritional data. Significant values (p < 0.05) for the Wilks’ lambda indicated in bold. 130 131 CHAPTER 4 Animal–mediated seed dispersal and the demo-genetic configuration across plant colonization gradients Jorge Isla, Miguel Jácome–Flores, Juan Miguel Arroyo, Cristina Rigueiro, Pedro Jordano & Cristina García (2024) In 2nd review in Journal of Ecology 132 133 ABSTRACT Ecologists have long recognized that seed dispersal mutualisms trigger natural regeneration and range expansion of animal-dispersed trees. Yet we lack empirical studies addressing whether frugivore activity influences founder effects, which reduce genetic diversity at the colonization front of expanding populations. Here we evaluate the contribution, from both demographic and genetic perspectives, of animal frugivores dispersing seeds across an expansion gradient. We used DNA barcoding for frugivores identification and highly polymorphic genetic markers (SSRs) for maternal analysis of juniper seeds to investigate how (1) stand maturity, (2) microhabitat types, and (3) foraging patterns shape the distribution of the maternal progenies along this gradient. We expect both reduced seed rain density and low numbers of source trees contributing to the seed rain at the colonization front, with limited availability of local fruiting trees. We also anticipated that large-sized frugivore species would promote maternally rich seed rain due to their ability to mix seed progenies during digestive processing and move further distances across the landscape. Contrary to our expectations, we found that all identified frugivores produced dense and genetically diverse seed rains across the expansion gradient, even at the colonization front, characterized by scarce fruiting trees. Our findings shed light on the fundamental and applied implications of plant-frugivore interactions in shaping highly diverse secondgrowth forests. These results emphasize the necessity of preserving plant-animal mutualistic interactions to ensure the persistence and expansion of natural tree populations, particularly in formerly fragmented landscapes. KEYWORDS: second–growth expanding forest, seed dispersal, seed rain, plant– animal mutualisms, maternal progenies, microhabitat, stand maturity, Mediterranean forests, frugivore’s forging patterns, natural regeneration. 134 Introduction Deforestation has led to a global reduction of forest cover over decades, but recent environmental policies and widespread rural land abandonment now foster forest regeneration and expansion worldwide (Chazdon et al. 2020). As secondgrowth forests expand, our dependence on them to support our wellbeing increases, but regrowth forests have long been considered less capable of providing ecosystem services comparable to those obtained from old-growth forests (Chazdon 2014). Yet, recent studies reveal that naturally regenerating forests on formerly tree-depleted rural lands exhibit comparable carbon uptake and host similar levels of plant and insect diversity (Espelta et al. 2020, Poorter et al. 1996). These results encourage further empirical studies to understand the regeneration dynamics underlying fast-paced forest natural expansion and, more importantly, to identify the factors that shape diverse, functional, and resilient future forests. Natural regeneration and the expansion of forest fragments commences with the formation of the seed rain, i.e., the spatial distribution throughout the landscape of all seeds dispersed in each reproductive event (Harper 1977, Figure 4.1). Seed rains, thus, serve as the starting template for plant regeneration in mature stands and for plant expansion in colonization front stands (Jordano & Godoy 2002). Frugivores provide dispersal services to animal-dispersed plants – including 70-94% of woody species in tropical forests and 50-75% in local Mediterranean habitats (Jordano 2014) – in exchange for food resources. Frugivore foraging patterns form strongly aggregated seed rain patterns, where seeds are disproportionately dispersed to few highly preferred deposition sites (see Figure 4.1 for further discussion on seed rain patterns, foraging behavior, and microhabitat preferences). The same foraging patterns also determine the distribution of maternal progenies within and across deposition sites, i.e., whether deposition sites receive seeds from one/few fruiting trees or, instead, they accumulate seeds from many different mother trees (Figure 4.1) (Karubian et al. 2010, García & Grivet 2011, Scofield et al. 2012, Carvalho et al. 2020). In this way, the demographic and genetic consequences of animal-mediated seed dispersal are 135 tightly intertwined. For example, Karubian et al. (2010) found that the foraging behavior of a long-wattled umbrellabird clumped maternal progenies of a palm tree at leks where males spend up to 80% of their time during their breeding season (see also Wenny & Levey 1998). Likewise, García et al. (2009) demonstrated that vegetation cover determined frugivores’ site preferences, with more maternal progenies found in covered microhabitats compared to open ones, a pattern corroborated in other studies (e.g., Lavabre et al. 2016). The spatial distribution of the maternal progenies across deposition sites is expected to influence the establishment of dispersed seeds and early stages of seedling recruitment (García et al. 2009, Grivet et al. 2005), but this aspect remains unexplored across population expansion gradients. For example, colonization fronts often exhibit sparse fruiting tree populations, potentially diminishing maternal source tree diversity in the seed rain while aggregating maternal progenies at a limited number of suitable deposition sites. 136 Seed shadows FORAGING PATTERNS SEED RAIN & SEED SHADOWS COMPONTENTS OF MATERNAL GENETIC CORRELATIONS Seed rain Maternal Richness (R) Maternal Redundancy (uPMI) Maternal Relatedness (r) A B C Overlapping seed shadows Non-overlapping seed shadow n = 8 R= 5n = 8 R= 2 n = 8 uPMI = 0.53 n = 8 uPMI = 0.25 n = 8 r= 0.70 n = 8 r= 0.10 Generalized foraging pattern Territorial foraging pattern 137 Figure 4.1. Distribution of maternal progenies in frugivore-generated seed rains. Frugivores transition from fruiting trees, where they forage on mature fruits, to deposition sites, where they rest, perch, or nest. In this process, they disperse seeds into deposition sites at preferred microhabitats (e.g. Jordano & Schupp 2000). Frugivores exhibit different foraging patterns defined by the combination of microhabitat preferences and foraging behavior that encompasses the sequence of movements among fruiting trees and deposition sites. (A) Territorial species typically display limited mobility, favoring specific preferred sites where dispersed seeds are clumped (A, orange line). In contrast, other species or individuals exhibit broader foraging behavior, traversing various fruiting-trees and widely distributed microhabitats throughout the landscape (A, blue line). (B) Frugivore foraging patterns imprint the seed rain, i.e., the spatial distribution of all seeds dispersed throughout the landscape from all fruiting trees that contribute seeds. The seed rain functions as the spatial template that shapes subsequent recruitment patterns of maternal progenies across the landscape, quantified as the number of seeds in a deposition site (expressed as seed density, no. seeds /m2). When the fruiting tree that produced dispersed seeds is known, we can depict the distribution of the maternal progenies throughout the landscape. This will show the extent of the seed shadows, which are the spatial distribution of all dispersed seeds from one specific mother tree (coloured shadows in B) (Janzen 1970). Overlapping seed shadows arise when different fruiting trees share deposition sites where they disperse progeny (B, green and orange seed shadows), while non-overlapping seed shadows occur when fruiting trees do not share deposition sites (B, blue seed shadow). (C) The maternal composition of the seed rain (distribution of progenies) is described using maternal genetic correlations at a given study unit that might include scats (C, brown ellipses), deposition sites (C, black polygons), microhabitats, or a nested combination of these study units (García et al. 2009). Indices describing the fine-scale genetic structure of these contributions to the seed rain include: Maternal Richness (R) refers to the number of distinct contributing source trees at the study unit (e.g. deposition site). Maternal Redundancy (uPMI) represents the maternal redundancy at the study unit and it is measured as the probability that two random seeds from the same deposition site come from the same tree (sensu Grivet et al. 2005). Maternal Relatedness (r) is the mean genetic relatedness (allele sharing probability) among the trees contributing seeds to a deposition site. Maternal relatedness increases when seeds come from nearby trees (seeds colored in green tones) and decreases when seeds come from distant trees because adult trees show strong spatial genetic correlations (Vekemans & Hardy 2004). 144 Figure 4.2. Seed rain of J. phoenicea across three stands in a population expansion gradient in Doñana National Park. A) Circles represent the locations of deposition sites with the size of the circle proportional to the estimated seed rain density and colored according to the microhabitat type. Black-ringed circles represent focal deposition sites in which contributing frugivores and dispersed progenies were identified. B) Variation in frugivore-mediated seed rain density in deposition sites among microhabitats in the three study stands (boxplots showing median, (50%) quartiles, range, and outliers). 145 Frugivore contribution to seed rain patterns across a heterogeneous expansion gradient We identified ten species of frugivores dispersing J. phoenicea seeds through DNA-barcoding of samples collected from focal deposition sites (N deposition sites= 105, N feces = 518, N seeds = 1000, Table 4.1, Figure 4.3). Turdus philomelos was the species that contributed most to the seed rain, accounting for 68.7% of dispersed seeds, followed by Erithacus rubecula, Turdus merula, and Vulpes vulpes, that collectively accounted for 28.4% of dispersed seeds (Table 4.1). Ancillary dispersers such as other Turdus species and non-legitimate dispersers (Chloris chloris and Cervus elaphus) contributed minimally to the seed rain (< 3%). Overall, frugivores contributed unevenly to stands (X2 = 95.3, d.f = 18, p < 0.001) and microhabitats (X2 = 365.3, d.f = 36, p < 0.001), (Figure 4.3ab). The main dispersers, T. philomelos and E. rubecula, contributed evenly to the seed rain across the three stands, while ancillary-species dispersed seeds mainly in mature and/or intermediate stands. All frugivores showed strong preferences for P. pinea and J. phoenicea microhabitats, where they consistently dispersed at least half of the seeds they mobilized, while only the main dispersers contributed to open area sites (Figure 4.3b). As a result, we found a wider range of frugivores contributing to P. pinea and J. phoenicea microhabitats, attracting the highest number of dispersed seeds, compared to microhabitats with a low number of dispersed seeds (open areas and other fleshy-fruited shrubs), (Figure 4.3c). The number of species contributing per deposition site ranged from 1, mostly in open area sites, to 5 under P. pinea. Turdus philomelos dispersed seeds to most deposition sites (Figure 4.3c, Figure C.1). 146 Table 4.1. Quantitative characterization of the seed rain and seed rain maternal composition sampled for J. phoenicea. Descriptors are indicated by frugivore species. Species listed in dark color are core species (dispersing > 97% of seeds) and species listed in grey are ancillary species (dispersing < 3% of seeds). 147 Figure 4.3. Relative contributions of the ten frugivore species to the seed rain of J. phoenicea at different spatial scales; (A) among stands, (B) among microhabitats and (C) among deposition sites. The relative contribution of frugivores per microhabitat is weighted by the number of deposition sites sampled. The species names of the 6 ancillary dispersers (that disperse <3% of the seeds) are gray colored. The bar Microhabitat cover represents the relative cover (%) of the microhabitats in the study area. Each column in C represents a deposition site, and the colors indicate the percentage of seeds found that have been dispersed by each frugivore species. The bottom, horizontal, bar represents the contribution of each frugivore to the overall seed rain. 148 Diversity and relatedness of maternal progenies in the seed rain Among 1000 dispersed seeds sampled, we identified 605 distinct maternal trees. Most maternal progenies arrive only at one stand, but two progenies reached both intermediate and colonization stands, proving the existence of LDD events (Figure C.2). These documented LDD events included at least two seeds (0.2% of the dispersal events assessed by seed genotyping) mediated by T. philomelos, resulting in dispersal distances ranging from 415-830 m. Within the three stands, deposition sites collected both seeds from unique progenies (only found at one deposition site) and progenies shared among deposition sites (mother trees that reached their progeny to multiple deposition sites), indicating overlapping of seed shadows (Figure C.2). Accumulation curves show that maternal richness increases with increasing sampled seeds across stands, microhabitats, and frugivore species, indicating that dispersed seeds come from a wide range of maternal trees (Figure C.3). Most maternal progenies were dispersed by T. philomelos (66.6%) followed by V. vulpes, E. rubecula, T. merula and ancillary-dispersers (12.2, 11.7, 9.7, and 4% respectively), reflecting the highly variable interaction frequencies of different frugivore species (Figure C.4). The maternal relatedness among seeds collected within the same scat was high compared to relatedness values of seeds sampled within and among deposition sites (Table 4.1). This result suggests that seeds carried within the same feces tend to come from the same or nearby source trees and that seeds collected in the same deposition site combine different maternal progenies. The components of maternal genetic correlations of the seed rain differed significantly among microhabitats (Figure 4.4, Table C.2), but not among stands (Figure C.5, Table C.2). Maternal richness differed significantly among microhabitats (F = 33.8, d.f = 4, p < 0.001) with P. pinea and J. phoenicea deposition sites showing the highest values. Mean relatedness and maternal redundancy (uPMI) also differed among microhabitats (F = 19.3, d.f = 4, p < 0.001; F = 6.1, d.f = 4, p < 0.001) with open area deposition sites showing increased values for both. This result suggests that 149 deposition sites in open areas, which receive a low number of seeds, tend to attract seeds from the same or closely related source trees. In contrast, all other deposition sites receive an even number of seeds from various poorly related source trees. CCA confirmed that the distinctive contribution (number of seeds) of this complex assemblage of frugivores explained the composition of maternal progenies in the seed rain. Up to 41% of the variation in diversity and relatedness of maternal progenies was explained by the differential contribution of frugivores (Wilk’s λ = 0.07, F = 23.3, d.f = 15,215, p < 0.001). The first canonical variable explained 99.5% of the total variance shared between the two original variable sets with a canonical correlation of R2 = 0.93 (F = 23.3, p< 0.001). Maternal richness was negatively correlated to this canonical variable (-0.99), while relatedness and redundancy were positively correlated (0.48 and 0.34 respectively) (Figure 4.5). The contributions of all frugivores were negatively correlated to this canonical variable with T. philomelos showing the highest correlation (-0.79). Although not significant, the second canonical variable was characterized by variations in the structure of relatedness among the deposition sites (Figure 4.5). These results confirm that species contributing with many seeds tend to increase the number of source trees and decrease relatedness and redundancy among progenies within deposition sites, resulting in extremely diverse seed shadows in terms of contributing maternal trees. The distribution of focal sites in the multidimensional CCA space, sorted by microhabitat type (Figure 4.5), but not by stand maturity (Figure C.6), shows that J. phoenicea and P. pinea sites tend to correlate with areas highly contributed by frugivores, where high levels of maternal diversity -and low levels of uPMI- were found. Deposition sites located under fleshy and non-fleshy fruited shrubs were scattered throughout the multidimensional space. On the contrary, most of the deposition sites in open areas were located in the space defined by low contributions of frugivores and reduced maternal richness, with increased values of relatedness and uPMI (Figure 4.5). 150 Figure 4.4. Maternal composition of the seed rain (Maternal richness, mean genetic relatedness and maternal redundancy uPMI) in deposition sites among microhabitats. Note that uPMI was log-transformed to improve visualization. Different letters indicate significant differences among microhabitat types. Boxplot showing median, (50%) quartiles, range, and outliers. Figure 4.5. Canonical correspondence analysis (CCA) biplot, showing the strength of association of seed contribution per frugivore species (arrows) with values of Maternal richness, relatedness and redundancy (uPMI) in deposition sites (dashed arrows) according to the length of the arrows. Animal species codes: Tphi = Turdus philomelos, Eru. = E. rubecula, Tmer = T. merula, Vvul = V. vulpes, An.sp. = Ancillary-species (C. chloris, C. elaphus, S. atricapilla, T. torquatus, T. iliacus and C. cooki). 151 Discussion Our results show that frugivores generate diverse seed rains regarding the maternal provenance of seeds across a colonization gradient. This is unexpected, according to the long-standing paradigm holding that founder effects, coupled with fast-paced population growth, diminish genetic diversity in the colonization front of expansion (Waters et al. 2013). Our methodological approach allowed us to connect the distribution of maternal progenies in the seed rain to the contributing frugivore species. This establishes a new perspective in the study of dispersal mutualisms, showing that: (1) a diverse assemblage of frugivores, including core and ancillary species, generates maternal-rich seed rains across an expansion gradient, even at the colonization front; and (2) contrary to our expectations, stand maturity does not impact the distribution of maternal progenies in the seed rain. Instead, the strong preferences of most active frugivores for certain microhabitats explains the consistent seed rain patterns regardless of stand maturity. We discuss the fundamental and applied outcomes of our main findings and suggest future research lines to advance our knowledge on the impact of plant-animal interactions underlying the regeneration dynamics and expansion of plant populations. Our study reveals that the accumulation of seeds and maternal progenies in the seed rain remains consistent across stands, irrespective of their maturity stage. We anticipated limited maternal contributions in areas with few local mother trees. However, contrary to our expectations, we observed a diverse array of maternal trees contributing progenies in the dense seed rain of the colonization front. This suggests thorough cone removal from most local trees, as well as not so infrequent long-distance seed dispersal. Our sampling provides evidence for LDD events (seed dispersal to different stands) for at least 0.2% of the 1000 seed dispersal events documented by endocarp genotyping. Most likely this represents an underestimate of the actual percentage of LDD events, given the amount of unique progeny arriving at each stand (Figure C.2). Yet, a large sample of adult trees must be genotyped to estimate LDD events frequency with more precision. 152 The recorded LDD events represent seed movements in the distance range of 415- 830 m (see Jordano 2017 for a discussion of within- and among-stand LDD events) and document the existence of among-stand seed dispersal. This includes seeds from external fruiting trees brought by frugivores in recently established stands (Hamrick & Trapnell 2011, Karubian et al. 2010, Lavabre et al. 2016). These results align with previous findings in the study area which noted a high number of visits by frugivores to fruiting trees at the colonization front (Isla et al. 2023). This process, even with relatively reduced (but not very rare) LDD events (Bialocyt et al. 2006), might contribute to maintaining moderate to high genetic diversity throughout the population expansion gradient (Browne & Karubian 2016, Scofield et al. 2012) mainly through the “reshuffling” and embolism effect involved in the LDD events (Bialocyt et al. 2006). Our study shows that on top of these non-trivial frequencies of LDD events, animal frugivores also contribute to the seed rain of very high diversity of provenances, even at the colonization front. Overall, these results highlight the potential of frugivores in mitigating the anticipated genetic erosion at the colonization front, which often arises due to founder effects followed by rapid population growth (Waters et al. 2013). This aligns with predictions from range expansion/colonization models where stratified dispersal patterns – combining both local in situ dispersal with LDD events – are expected to yield favorable conditions preserving genetic diversity by swamping the effects of drift and the homogenizing effect of surfing (Bialocyt et al. 2006, Excoffier et al. 2009, Le Corre & Kremer 1998). Our results show that here, the mutualistic interaction with a diversified frugivore assemblage provides these conditions and that a highly structured seed rain is the primarily outcome of such interaction. This baseline study presents a new perspective on the study of dispersal mutualisms and may prompt further research to explore whether the potential of this interaction for promoting diverse seed rains is sustained in different ecological systems during natural expansion. We found a frugivore assemblage consisting of birds and mammals, resembling those previously documented in juniper populations (Escribano-Ávila et al. 2014, Jordano 1993, Livingston 1972, Santos et al. 1999). This assemblage 153 consisted of a core of four primary frugivores, predominantly T. philomelos that contributed 68.7% of seeds to the seed rain, along with six ancillary species. The key role of thrushes as seed dispersers is widely acknowledged across various animal-dispersed trees (Escribano-Ávila et al. 2014, Jordano 1993, Rumeu et al. 2009, Snow & Snow 1988). Furthermore, our data reveals that T. philomelos primarily disperse seeds of J. phoenicea, revealing a symmetric interaction between this thrush and juniper trees (Bascompte et al. 2006). This symmetry is likely driven by the phenological alignment of both species in the area, as J. phoenicea is the dominant fruiting species during the wintering period of T. philomelos. This mutual reliance signifies a reciprocal beneficial relationship (Quintero et al. 2023), that presumably sustains migrant populations of keystone frugivore species in the area, while concurrently enhancing the regeneration and expansion dynamics of the locally abundant juniper tree. Frugivore preferences for specific microhabitats and their foraging behavior are the two key elements shaping the density and maternal composition patterns in seed rains (García et al. 2009, Jordano & Schupp 2000, Tewksbury et al. 2002). The initial determinant entails frugivore preference for distinct microhabitats. These preferred microhabitats are characterized by dense vegetation that facilitates feeding, perching, and predator avoidance sites (e.g. P. pinea, J. phoenicea and fleshy-fruited shrubs in our study), contrasting with frugivore elusiveness to low or non-vegetated sites (non-fleshy fruited shrubs and open areas in our study) (García et al. 2009, Jordano & Schupp 2000). These strong microhabitat preferences infuse seed rain patterns, as highly frequented sites amass a broad range of maternal progenies, reflecting the foraging preferences of the most active frugivores (García et al. 2009, Karubian et al. 2010). Corroborating this, we found high seed density and maternal richness in P. pinea and J. phoenicea deposition sites. Our findings reveal that fewer than 20% of seeds beneath juniper trees are maternal half-sibs (i.e., same mother tree, Figure C.2) contrasting with previous studies where 70% of the animal-dispersed seeds were deposited under the mother tree (Jordano & Godoy 2002). These findings reveal frugivores’ transitory foraging behaviors among multiple fruiting trees before dispersing seeds