Interactions, Heterogeneity, and the Determinants of Multispecies Coexistence

Universidad de Sevilla Departamento de Biolog´ ıa Vegetal y Ecolog´ ıa DOCTORADO EN BIOLOG´ IA INTEGRADA Interactions, Heterogeneity, and the Determinants of Multispecies Coexistence Rodrigo Rodr´ ıguez Granjel Advisors: Dr. Oscar Godoy del Olmo Dr. Juan Arroyo Mar´ ın Tutor: Dr. Juan Arroyo Mar´ ın February, 2023 i Rodrigo R. Granjel was supported by a ”Formaci ´ on de Profesorado Universitario” (FPU) contract (2017—2022) from the former Spanish Ministry of Education, Culture and Sports and by a Fulbright fellowship (2022—2023) from the Spanish Fulbright Commission and the Spanish Ministry of Foreign Affairs, sponsored by the Junta de Andaluc ´ ıa. The author developed the current doctoral thesis as part of the Ecology, Evolution, and Conservation of Mediterranean Plants research group (EVOCA), the Mediterranean Ecological Studies and Synthesis laboratory (MESSY Lab), and the Radical Community Ecology group (RCE). Specific project funding information is provided at the end of each corresponding chapter. The different studies were carried out in, or have greatly benefited from, the following research institutions: •Universidad de Sevilla (US; Sevilla, Spain). FPU fellowship (2017—2022). • Instituto de Recursos Naturales y Agrobiolog ´ ıa de Sevilla (IRNAS-CSIC; Sevilla, Spain). Training agreement US-IRNAS (2017—2018). • Institute of Plant Sciences, University of Bern (IPS-UNIBE; Bern, Switzerland). Mixed funds (2017, 2018). Swiss-European Mobility Programme, StaffMobility Training (2019). • Universidade Estadual Paulista (UNESP; Rio Claro, Brazil). Santander Iberoam´ erica Investigaci´ on fellowship (2019—2020). • Misi ´ on Biol ´ oxica de Galicia (MBG-CSIC; Pontevedra, Spain). FPU mobility grant (2021). • North Carolina State University (NCSU; Raleigh, North Carolina, USA). Fulbright predoctoral fellowship (2022—2023). ii iii Statements for each chapter regarding associated publications, code and data availability, additional funding, and author contributions, are provided at the end of each corresponding section. Credits for the illustrations on the cover of each chapter: • Introduction — Monet, Claude. Banks of the Seine. 1880. National Gallery of Art, Washington, DC, USA • Chapter 1 — Sisley, Alfred. Meadow. 1875. National Gallery of Art, Washington, DC, USA • Chapter 2 — van Gogh, Vincent. Plain near Auvers. 1890. Neue Pinakothek, Munich, Germany • Chapter 3 — van Gogh, Vincent. Olive trees. 1889. The Nelson-Atkins Museum of Art, Kansas City, MO, USA • Chapter 4 — Haes, Carlos de. Marismas (Holanda). ca. 1884. Museo del Prado, Madrid, Spain • Discussion — Corot, Camille. The Ferryman. 1865. The Metropolitan Museum of Art, New York, NY, USA iv v Abstract No living being is an island. Every creature needs resources from its environment to survive and thrive throughout its life. But in a finite world with over 8.7 million species, acquiring such resources is an inherently competitive endeavour. This has resulted in an extraordinary diversity of interactions, with varying strengths that have either negative or positive consequences for species’ performance. Over the last century, ecologists have sought to understand how the balance of these interactions affects a species’ ability to persist in a community. Initially, research focused on studying interactions between pairs of species, finding that coexistence is possible when stabilising niche differences, which occur when intraspecific competition exceeds interspecific competition, outweigh fitness differences that favour one competitor over others. However, species typically engage in interactions that go beyond species pairs, making this pairwise framework unsuitable to explore species persistence in complex communities. Besides, one main characteristic of ecological interactions is that they are not set in stone—on the contrary, they vary through space and time in response to modifications in the surrounding conditions. Therefore, new theoretical developments are providing us with tools to comprehend how interactions shape biodiversity in complex, multispecies communities with environmental heterogeneity, which is of paramount importance in the changing world that we inhabit. In this thesis, I combine the most recent advances in ecological theory with thorough empirical data to explore how different sources of environmental heterogeneity alter biotic interactions between plant species, and in turn how the structure of such interactions shapes the determinants of multispecies coexistence. In the first three chapters, I use a structural stability approach that allows me to investigate the mechanisms of coexistence beyond pairwise combinations of species. In chapter one, I examine how changes in resources (nitrogen availability) and natural enemies (foliar pathogens) modify the mechanisms of plant diversity and composition. To do that, I quantify the intrinsic growth rates and interaction coefficients between eight common Central European perennial plants in an experiment factorially modifying the nitrogen availability and foliar fungal pathogens of a Swiss perennial grassland. I find that both nitrogen addition and pathogen suppression decrease structural fitness differences that drive competitive dominance but, surprisingly, they also promote niche differences that stabilise the dynamics of interacting species. Interestingly, all effects of resources and enemies on the mechanisms of plant coexistence are dependent on the number of interacting species. In chapter two, I use data from 8 years and 150 German grasslands to investigate the effects of land use intensification on the growth rates and interactions of 50 perennial plant species. The data shows that increasing land use causes species loss by reducing structural niche differences in a non-linear way, rather than by increasing differences in fitness. However, I also find that niche differences play a role in maintaining coexistence among the species that persist at high land use intensities. xii CONTENTS 3.2.1 Studysystem........................... 70 3.2.2 Grasshoppers and experimental design . . . . . . . . . . . . 70 3.2.3 Experimental procedure . . . . . . . . . . . . . . . . . . . . 70 3.2.4 Interaction coefficients and structural stability . . . . . . . 71 3.2.5 Contribution of each interaction layer . . . . . . . . . . . . 72 3.2.6 Modulemetrics ......................... 74 3.3 Results .................................. 76 3.4 Discussion ................................ 81 4 Precipitation, traits, and niche variability 87 4.1 Introduction ............................... 91 4.2 Materials and methods . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.2.1 Study system and fieldwork . . . . . . . . . . . . . . . . . . 93 4.2.2 Estimating species interactions . . . . . . . . . . . . . . . . 94 4.2.3 Calculation of niche and fitness . . . . . . . . . . . . . . . . 94 4.2.4 Functional traits and climatic variables . . . . . . . . . . . . 96 4.2.5 Statistical analyses . . . . . . . . . . . . . . . . . . . . . . . . 96 4.3 Results .................................. 98 4.4 Discussion ................................ 103 Discussion and conclusions 107 Appendix 1 — Chapter 1 121 Appendix 2 — Chapter 2 159 Appendix 3 — Chapter 3 169 Appendix 4 — Chapter 4 181 Bibliography 189 List of Figures 0.1 Structural stability approach . . . . . . . . . . . . . . . . . . . . . . 11 1.1 Environmental changes and structural coexistence . . . . . . . . . 23 1.2 Experimental design in PaNDiv . . . . . . . . . . . . . . . . . . . . 27 1.3 Example of structural coexistence . . . . . . . . . . . . . . . . . . . 32 1.4 Nitrogen, pathogens, and niche and fitness differences . . . . . . . 33 1.5 Indirect effects of nitrogen and pathogens . . . . . . . . . . . . . . 34 1.6 Nitrogen, pathogens, and feasibility . . . . . . . . . . . . . . . . . . 35 1.7 Nitrogen, pathogens, and community composition . . . . . . . . . 36 2.1 LUI and structural niche and fitness differences . . . . . . . . . . . 56 2.2 Indirect effectsofLUI.......................... 59 2.3 Coexistence observations and predictions with LUI . . . . . . . . . 60 3.1 Scenarios of network complexity . . . . . . . . . . . . . . . . . . . 73 3.2 Empirical network and interaction distributions . . . . . . . . . . 76 3.3 Contribution of interactions to coexistence . . . . . . . . . . . . . . 77 3.4 MetricsandSND ............................ 78 3.5 MetricsandSFD............................. 79 4.1 Niche and fitness, ellipses per year . . . . . . . . . . . . . . . . . . 98 4.2 Species niche and fitness ranking per year . . . . . . . . . . . . . . 99 4.3 Species’ fitness and precipitation . . . . . . . . . . . . . . . . . . . 101 4.4 Species’ niche and precipitation by season . . . . . . . . . . . . . . 102 S1.1 Biomass-cover relationships . . . . . . . . . . . . . . . . . . . . . . 122 S1.2 Biomass-cover relationships . . . . . . . . . . . . . . . . . . . . . . 123 S1.3 Structural niche differences per year . . . . . . . . . . . . . . . . . 124 S1.4 Structural niche differences 95% confidence . . . . . . . . . . . . . 125 S1.5 Structural fitness differences per year . . . . . . . . . . . . . . . . . 126 S1.6 Structural fitness differences 95% confidence . . . . . . . . . . . . 127 S1.7 Indirect effectsperyear......................... 128 S1.8 Indirect effects 95% confidence . . . . . . . . . . . . . . . . . . . . 129 S1.9 Feasibility per year . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 S1.10 Feasibility 95% condifence . . . . . . . . . . . . . . . . . . . . . . . 131 S2.1 LUIvaluesperyear ........................... 160 S2.2 LUI and species interactions . . . . . . . . . . . . . . . . . . . . . . 161 S2.3 Quantile regressions SND, SFD and LUI . . . . . . . . . . . . . . . 162 S2.4 Intrinsic growth rate and interspecific interactions in LUI . . . . . 163 xiii xiv LIST OF FIGURES S2.5 SND and SFD across LUI, lower 95% CI . . . . . . . . . . . . . . . 164 S2.6 SND and SFD across LUI, upper 95% CI . . . . . . . . . . . . . . . 165 S3.1 Experimental design . . . . . . . . . . . . . . . . . . . . . . . . . . 171 S3.2 Cross validation error and lambda . . . . . . . . . . . . . . . . . . 172 S3.3 Contribution of interactions to coexistence . . . . . . . . . . . . . . 173 S3.4 Metrics and SND, 4 species . . . . . . . . . . . . . . . . . . . . . . . 174 S3.5 Metrics and SND, 4 species . . . . . . . . . . . . . . . . . . . . . . . 175 S4.1 Species’ fitness and precipitation by season . . . . . . . . . . . . . 182 List of Tables 4.1 PERMANOVA functional traits . . . . . . . . . . . . . . . . . . . . 100 S1.1 List of plant species in PaNDiv . . . . . . . . . . . . . . . . . . . . 132 S1.2 AIC and model selection in PaNDiv . . . . . . . . . . . . . . . . . . 133 S1.3 Control αand λcoefficients ...................... 134 S1.4 Nitrogen αand λcoefficients ..................... 135 S1.5 Fungicide αand λcoefficients..................... 136 S1.6 Combined αand λcoefficients .................... 137 S1.7 GLMsummaryofSND......................... 138 S1.8 GLMsummaryofSFD ......................... 139 S1.9 GLM summary of SND in 2017 . . . . . . . . . . . . . . . . . . . . 140 S1.10 GLM summary of SND in 2018 . . . . . . . . . . . . . . . . . . . . 141 S1.11 GLM summary of SFD in 2017 . . . . . . . . . . . . . . . . . . . . . 142 S1.12 GLM summary of SFD in 2018 . . . . . . . . . . . . . . . . . . . . . 143 S1.13 GLM summary of SND, lower 95% . . . . . . . . . . . . . . . . . . 144 S1.14 GLM summary of SFD, lower 95% . . . . . . . . . . . . . . . . . . 145 S1.15 GLM summary of SND, upper 95% . . . . . . . . . . . . . . . . . . 146 S1.16 GLM summary of SFD, upper 95% . . . . . . . . . . . . . . . . . . 147 S1.17 GLM summary of differential ..................... 148 S1.18 GLM summary of differential in 2017 . . . . . . . . . . . . . . . . 149 S1.19 GLM summary of differential in 2018 . . . . . . . . . . . . . . . . 150 S1.20 GLM summary of differential, lower 95% . . . . . . . . . . . . . . 151 S1.21 GLM summary of differential, upper 95% . . . . . . . . . . . . . . 152 S1.22 GLM summary of feasibility . . . . . . . . . . . . . . . . . . . . . . 153 S1.23 GLM summary of feasibility in 2017 . . . . . . . . . . . . . . . . . 154 S1.24 GLM summary of feasibility in 2018 . . . . . . . . . . . . . . . . . 155 S1.25 GLM summary of feasibility, lower 95% . . . . . . . . . . . . . . . 156 S1.26 GLM summary of feasibility, upper 95% . . . . . . . . . . . . . . . 157 S2.1 LUI quantile regression statistics . . . . . . . . . . . . . . . . . . . 166 S2.2 LUI quantile regression statistics, lower 95% CI . . . . . . . . . . . 167 S2.3 LUI quantile regression statistics, upper 95% CI . . . . . . . . . . 168 S3.1 Plant species in the grasshoppers experiment . . . . . . . . . . . . 176 S3.2 GLM output for SND and module metrics . . . . . . . . . . . . . . 177 S3.3 GLM output for SFD and module metrics . . . . . . . . . . . . . . 178 S3.4 GLM for SND and module metrics, 4 species . . . . . . . . . . . . 179 S3.5 GLM for SFD and module metrics, 4 species . . . . . . . . . . . . . 180 S4.1 Plant species list from Caracoles . . . . . . . . . . . . . . . . . . . 183 S4.2 Pairwise t-test, fitness and year . . . . . . . . . . . . . . . . . . . . 184 xv xvi LIST OF TABLES S4.3 Pairwise t-test, niche and year . . . . . . . . . . . . . . . . . . . . . 185 S4.4 GLM fitness and precipitation . . . . . . . . . . . . . . . . . . . . . 186 S4.5 GLM niche and precipitation . . . . . . . . . . . . . . . . . . . . . 187 Introduction 1 Introduction 2 Biodiversity maintenance: Solving the paradox of the plankton Planet Earth is home to over 8.7 million living species of bacteria, archaea, protists, fungi, plants, and animals (Mora et al., 2011), yet only 1.2 million of them have been catalogued to this date. Such diversity of life has fascinated humanity for centuries, but it was not until the advent of ecology as a distinct discipline that we began to uncover the mechanisms by which multiple species are able to coexist. Ecologists began by trying to explain the observed natural associations between species and environmental variation. One of the foundational concepts in community ecology is the species’ niche, first introduced by Joseph Grinnell (1917) as the range of environmental conditions under which a species could survive. But this concept, far from uniform, has a long history of developments and mutations. Soon after Grinell’s definition, Charles Elton (1927) independently defined the niche as the function or role that a species plays in a community, mainly through interactions with food and predators. Not much later, Georgy Gause (1934) proposed a connection between both abiotic and biotic-functional definitions of the niche. This inspired George E. Hutchinson (1957) to conceptualise the niche as a multi-dimensional space formed by all the environmental factors that a species endures, which is the base of our current understanding of this fluctuating, complex ecological concept. In parallel with the development of the niche concept, a competitive framework for community ecology was being gestated. Together with the formulation of Alfred Lotka (1925)’s and Vito Volterra (1926)’s famous mathematical models on competition between two different species (i.e., interspecific competition), came the prediction that two species competing for a single resource may not coexist. This prediction was later experimentally demonstrated by Gause (1934), who showed that two protists growing on a single resource in controlled conditions can achieve stable populations in isolation but one of them excludes each other in combination. This phenomenon of two competing species being unable to coexist on a single resource was subsequently known as Gause’s competitive exclusion principle. With both these concepts in mind, Robert H. MacArthur and Richard Levins (MacArthur and Levins, 1967) attempted to understand the coexistence of competing species by expanding on the Hutchinsonian niche and proposing the concept of limiting similarity, which can be defined as the minimal difference in species niche requirements that allows them to coexist; or else, the maximum niche overlap. By simplifying the multidimensional space and focusing on a few critical niche axes in which competition occurs (MacArthur 1969), this approach showed a way to link the niche overlap to the processes of competitive exclusion in a quantitative manner (MacArthur, 1972), in both deterministic and fluctuating environments (May and MacArthur, 1972; May, 1973). These ideas conceived species coexistence mainly through competition for limited resources: two species needed to display enough niche differences to reduce their interspecific competition, given that a high niche overlap would lead to the competitive exclusion of one of the species. Therefore, an ecological system could support as many species as there were available niches, i.e., n different factors or resources could accommodate n species. However, this led to a paradox, 3 Introduction Therefore, the 2-species structural fitness differences (SFD) can be defined as the deviation of r from rc , so the closer r is to the optimum set of intrinsic growth rates determined by the structure of species interactions, the lesser the SFD: SFD = arccos r·rc ||r||·||rc||!(0.7) This approach can be applied to all three pairs found in our three-species community. To simplify, the feasibility domains of the three pairs of species can be visualised in a unit simplex triangle (Fig. 0.1c), which is a generalisation with unit area to n-dimensions. In this triangle, each side represents the full domain of intrinsic growth rates between the corresponding vertex species in such a way that allows us to project the feasibility domain between the pair of species involved. The pedagogical advantage of this simplex representation is that, by projecting the uni-dimensional feasibility domain of each one of the three pairs (light green cones) into the two-dimensional inner parameter space of the triangle, one can easily observe the range of 3-species intrinsic growth rate vectors compatible with different feasible equilibrium points. Importantly, the overlap of the projections of the different feasibility domains signals the common set of intrinsic growth rates compatible with feasible combinations of species, and may include a common domain in which all three domains overlap called pairwise feasibility domain (darker shade of green). As an example, in Fig. 0.1c, the 3-species intrinsic growth rates vector, r (orange dot), lies within the overlapping region of the feasibility domains of pairs ik and jk . Extending the feasibility condition in eqn. 0.4 and as seen in Fig. 0.1b, this means that there are chances for the coexistence of all three species if they interact in ik and jk pairs, but the pair ij alone would not be able to achieve a feasible equilibrium provided the current intrinsic growth rates. But the main contribution of the structural stability approach is the ability to calculate SND and SFD for n-dimensional systems. Fig. 0.1d shows the full positive parameter space of intrinsic growth rates for our three-species community. Similarly to the 2-species case, the structure of the matrix of species interactions, α , determines a range of feasible solutions that can accommodate the differences in intrinsic growth rates, r . Calculations shown in eqn. 0.5 through eqn. 0.7 are for the specific case of a 2-species community; the structural coexistence mechanisms generalised for nnumber of species can be calculated as follows: SND = |det(α)| n √π/2Z···ZRn ≥0 e−xTαTαxdx. (0.8) rc=1 nvi ∥vi∥+vj vj+···+vn ∥vn∥.(0.9) 10 Figure 0.1: Graphical representation of the structural stability approach to multispecies coexistence. a, three-species ( i , j , and k ) interaction network, with node size representing the strength of intrinsic growth rates ( ri , rj , and rk ) and link width representing the strength of intraspecific ( αii , αjj , and αkk ) and interspecific ( αji , αki , αij , αkj , αik , and αjk ) interaction coefficients. (Continued.) 11 Introduction Figure 0.1:(Continued.) b, graphical description of the calculation of structural niche differences (SND) and structural fitness differences (SFD) for the pair ij in panel a. The dark green slopes αjj/αij and αii/αji limit the feasibility domain, which corresponds to the SND or the normalised solid angle between both slopes (dark green curved line). The SFD can be calculated as the angle between the centroid ( rc ) or the feasibility domain and the vector formed by the intrinsic growth rates of both species ( r ). c, geometric projection on the unit simplex for the feasibility domains (light green triangles) between the three species in a(e.g., between i and k as pointed in the panel). The overlap between the three feasibility domains creates a region called pairwise feasibility domain, which indicate the range of intrinsic growth rates ( r ) that are compatible with the coexistence of the three pairs of species. In this example, r lies on the region where the feasibility domains of jk and ik overlap, determining the feasibility for these pairs but not for the pair ij .d, Expansion of the pairwise structural approach for multispecies communities. Instead of two axes, we now have a parameter space of n axes formed by the r of the n species involved (in this example, the three species in panel a). As seen in panel b, the structure of species interactions defines a multidimensional parameter space of intrinsic growth rates compatible with the feasibility of the n species, the multispecies feasibility domain (green hypervolume), the solid angle of which corresponds to the measurement of multispecies SND. SFD corresponds to the angle between the vector of intrinsic growth rates, r , and the centroid of the feasibility domain rc . e, geometric projection on the unit simplex for the three-species feasibility domain, with indication of the approximation of SND and SFD in this representation. In this example, the three species are not predicted to coexist as the SFD are larger than the SND ( r does not lie inside of the 3-species feasibility domain). Comparing the size and shape of the multispecies feasibility domain with the pairwise feasibility domain can be a powerful tool to understand the role of indirect interactions, which only emerge in multispecies assemblages, to promoting coexistence. Finally, fprovides a representation of the relationship between SND and SFD in a number of random 3-species communities, together with their feasibility predictions. The green triangle represents the region of SND and SFD values for which a pairwise combination is predicted feasible. Comparing the 3-species feasibility predictions with the 2-species feasibility predictions allows us to understand the importance of indirect interactions, as the limit that determines feasibility is not as clear for 3 as it is for 2 species. Panels b—eafter Saavedra et al. (2017). SFD = arccos r·rc ||r||·||rc||!.(0.10) Using the same simplex triangle as for the 2-species approach (Fig. 0.1c), we can now visualise the feasibility domain for the triplet (dark green; Fig. 0.1e). Note that, in this example, the multispecies feasibility domain differs with respect to the pairwise feasibility domain. The only reason this happens is the effect of indirect interactions that only emerge in multispecies systems, such as intransitivity (Saavedra et al., 2017; Godoy et al., 2017; Soliveres et al., 2015). The best example of an intransitive loop of indirect interactions is the well-known game of rock-paper-scissors, where all components lose and win against one other component of the network. Take only two of the three and one will exclude the other, but take the three of them together and they reach a feasible equilibrium. Therefore, through the structural stability approach it is possible to quantify the contribution of indirect effects to multispecies coexistence by comparing the pairwise and the multispecies feasibility domains through two different metrics. First, by comparing the relative size of the multispecies versus the pairwise feasibility domain, we are able to understand whether indirect interactions hinder 12 or help promote coexistence opportunities. This metric, named community-pair differential by Saavedra et al. (2017) (hereafter, differential), takes continuous values from -1 to 1, with negative values indicating less and positive values indicating more coexistence opportunities in multispecies assemblages compared to the pairwise approach. Second, by analysing the percentage of the multispecies feasibility domain that lies inside of the pairwise domain we can understand the contribution of pairwise mechanisms to explaining multispecies coexistence opportunities. This metric was named community-pair overlap (hereafter, overlap), and takes values from 0 (when there is no overlap between the two domains) to 1 (when the whole multispecies feasibility domain lies inside the pairwise domain). Note that, although both metrics are correlated (for example, the higher the differential, the lower the overlap will tent to be), they describe well differentiated mechanisms driven by indirect interactions. Finally, Fig. 0.1f shows feasibility predictions for a number of 3-species communities, with random α and r coefficients, in relation to their calculated multispecies SND and SFD. Additionally, the green triangle in the background of the panel indicates the feasibility condition for any 2-species community (i.e., any pair of species showing SND and SFD values within this area would be predicted feasible, and vice versa). By comparing how the 2-species and the 3-species behave to determine coexistence, we observe that, unlike the hard transition from feasibility to unfeasibility seen for species pairs, the 3-species communities show a transition region in which the feasibility prediction is uncertain and seemingly independent of the structural niche and fitness differences. This phenomenon is evidence of the importance of indirect interactions on determining the feasibility of multispecies communities. Further details on the structural stability approach can be found in Rohr et al. (2014), Saavedra et al. (2017), and Song et al. (2018). For a graphical illustration of this approach, I strongly recommend the interactive shiny app by Petry and Lepori (2022), available at the following web address as of February, 2023: https://ecodynamics.shinyapps.io/StructuralCoexistence. Importantly, the R code to calculate the multispecies mechanisms of species coexistence is available through Saavedra et al. (2017) and also through Granjel et al. (2022). To date, most advances around the structural coexistence framework remain theoretical or apply a combination of theory and demonstrations using existing empirical datasets (Cenci and Saavedra, 2018, 2019; Song and Saavedra, 2018; Song et al., 2020b; Saavedra et al., 2020; Medeiros et al., 2023). For example, in their original paper, Saavedra et al. (2017) used a dataset previously analysed under modern coexistence metrics (Godoy and Levine, 2014; Godoy et al., 2014; Kraft et al., 2015a) to show that coexistence predictions increase under a structural perspective. Only a handful of ad hoc empirical studies have been carried on. One of the first results was found by Petry et al. (2018), who showed that a competition-defence trade-off, where ants preferentially consume bigger, competitive dominant seeds in an annual plant community with stabilising niche differences, did not increase diversity but reduced coexistence by disadvantaging and even excluding the superior competitor through overconsumption, particularly in multispecies assemblages. Garc ´ ıa-Callejas et al. (2021) explored the role of spatial variability in the coexistence of an annual marshland, finding that 13 Introduction spatial heterogeneity plays little role in the accumulation of species richness with increasing spatial scale; instead, more interactions are realised and species coexist through direct rather than indirect effects. Working with microbes, Deng et al. (2021) found that non-resident species are likely to change the game rulebook of resident communities, establishing and suppressing resident species under changing environments. Also, Tabi et al. (2020) revealed that, at warmer temperatures, microbial communities show relatively more asymmetry in their feasibility domain and idiosyncrasy of responses compared with those communities at cooler temperatures. In a resolute multitrophic manipulation experiment, Bartomeus et al. (2021) showed that the persistence of plants increases when accounting for plant-pollinator interactions and that multitrophic structure manipulation results in unpredicted changes in persistence probabilities. Johnson et al. (2022), however, found that competition for pollinators destabilises interactions and plant coexistence, and competitive imbalances are expected to increase as pollinators are lost, hindering the coexistence of assemblages of 3, 4 and 5 species. Also for pollinators, recently Arroyo-Correa et al. (2023) found that communities with high plant species-level degree and low overlap among species are more feasible and that interindividual variation plays an important role in promoting feasibility by lessening the overlap in pollinator use among plant species, with communities composed of variable populations being 28% more feasible than communities comprising similar individuals. Brief statement on knowledge gaps and research opportunities Although hundreds of different theoretical and empirical studies emerged after modern coexistence theory was formalised with Chesson (2000), and therefore we have acquired a well-rooted understanding of how abiotic and biotic heterogeneity shapes the mechanisms of coexistence between pairs of species, the current state of the literature reveals clear gaps in our knowledge on how coexistence mechanisms operate in multispecies, more realistic species communities (Levine et al., 2017). Recent efforts have shed light on different key aspects of the ecology of natural communities, but we are lacking studies on the role of main drivers of diversity in ecological systems, such as fluctuations in resources and habitat availability, natural enemies, or temporal variability, in determining the fate of multiple interacting species. Overall, the latest theoretical developments in community ecology offer valuable tools to comprehend the mechanisms of species persistence and biodiversity maintenance in complex ecosystems, and enable us to connect these findings to the underlying abiotic and biotic heterogeneity in the face of global change. 14 Objectives The general goal of this thesis is to combine the most recent developments in ecological theory with thorough empirical data to begin to understand how species interactions shape the mechanisms of species coexistence in response to different sources of biotic and abiotic environmental heterogeneity. I address fundamental questions regarding the consequences for biodiversity of resource availability, natural enemies, land use intensity, different types of interactions, and temporal variability. Specifically, I propose the following objectives: • In Chapter 1, I aim to examine how changes in resources (nitrogen availability) and natural enemies (foliar pathogens) modify the mechanisms of plant diversity and composition. Specifically, I attempt to elucidate whether the mechanisms of multispecies plant coexistence differ in their response to resource availability and natural enemies from the mechanisms of coexistence between pairs of species. Also, I seek to measure the contribution of indirect interactions between plant species to explain the multispecies mechanisms of coexistence. • In Chapter 2, I seek to investigate the effects of land use intensification on the mechanisms of multispecies coexistence. In particular, I intend to understand whether the pattern of species loss that generally occurs following land intensification is explained by a reduction in structural niche differences that promote coexistence through self-limitation or by an increase in the structural fitness differences that favour one species over the others. • In Chapter 3, I propose to evaluate the role of different components of interaction networks in shaping structural niche and fitness differences. In detail, I aim to analyse the contribution of direct interactions between plants, higher-order interactions between plants, direct interactions between plants and herbivores such as grasshoppers, and higher-order interactions between plants and grasshoppers, to the multispecies mechanisms of plant coexistence. Moreover, I attempt to describe how different metrics of interaction configurations (sign, strength distribution, and intraspecific vs. interspecific effects) predict the observed structural niche and fitness differences between plant species. • In Chapter 4, I attempt to apply the latest developments in modern coexistence theory to characterise niche and fitness at the species level, rather than calculating differences between interacting species. With this in mind, I seek to describe the temporal variation in a plant’s niche and fitness and understand whether such variation can be explained by species-specific features (e.g., functional traits) or whether it responds better to environmental characteristics (e.g., rainfall). 15 Nitrogen enrichment and foliar fungal pathogens affect the mechanisms of multispecies plant coexistence 1 17 1. Nitrogen, pathogens, and multispecies coexistence 18 Abstract Changes in resources (e.g., nitrogen) and enemies (e.g., foliar pathogens) are key drivers of plant diversity and composition. However, their effects have not been connected to the niche and fitness differences that determine multispecies coexistence. Here, we combined a structural theoretical approach with a detailed grassland experiment factorially applying nitrogen addition and foliar fungal pathogen suppression to evaluate the joint effect of nitrogen and pathogens on niche and fitness differences, across a gradient from two to six interacting species. Nitrogen addition and pathogen suppression modified species interaction strengths and intrinsic growth rates, leading to reduced multispecies fitness differences. However, contrary to expected, we also observed that they promote stabilising niche differences. Although these modifications did not substantially alter species richness, they predicted major changes in community composition. Indirect interactions between species explained these community changes in smaller assemblages (three and four species) but lost importance in favour of direct pairwise interactions when more species were involved (five and six). Altogether, our work shows that explicitly considering the number of interacting species is critical for better understanding the direct and indirect processes by which nitrogen enrichment and pathogen communities shape coexistence in grasslands. Keywords Grasslands, nitrogen enrichment, foliar fungal pathogens, multispecies coexistence, structural approach, niche differences, fitness differences, indirect effects 19 1. Nitrogen, pathogens, and multispecies coexistence plant species in temperate grasslands and the selection involves a mix of grasses and herbs and fast- and slow-growing species. For each species, we measured the biomass of multiple individuals in different competition neighbourhoods to calculate the strength and sign of species interactions and how they varied across nitrogen and fungicide treatments. The natural variation in density, and the experimentally imposed relative frequency of these species across plots allowed us to separate and estimate with confidence intra- and interspecific competition as well as species’ intrinsic growth rates (see below). We then parameterised a model that describes the population dynamics of competing species, and their likelihood of coexisting, as a function of their network of interactions and their intrinsic growth rates in the absence of competition (see below, structural approach). With eight species, we need to estimate 64 pairwise interactions per treatment (8 intraspecific interaction coefficients and 56 interspecific) and a total of 256 pairwise interactions across the four treatments (i.e., control, nitrogen only, fungicide only, or both combined). To account for this high dimensionality, we proceeded as follows. Prior to the growing season of 2017 (February and early March), we marked focal individuals of each species in neighbourhoods varying in density and relative frequency of co-occurring species, using a bamboo stick together with a coloured plastic tag. We defined the neighbourhood for each focal individual as a 40cm diameter circle, marked with a PVC ring during sampling. We selected neighbourhoods in the monoculture plots to estimate intraspecific interactions and we did the same in the polyculture plots to estimate interspecific interactions (Fig. 1.2). We measured eight randomly selected neighbourhoods per pairwise interaction, across the four blocks in order to produce gradients of neighbour density (8 neighbourhoods/pairwise interaction × 256 pairwise interactions in the whole experimental design = 2,048 neighbourhoods). A minimum cover of at least 5% of either the same or another species was the criterion taken to measure a selected neighbourhood. Those neighbourhoods that did not meet this criterion (because they were dominated by some of the other 12 species in the experiment) were replaced by another one within the same plot. For example, a focal-neighbour combination could be T. officinale surrounded by at least 5% cover of A. odoratum. It is worth noting that, contrary to previous studies in which polycultures were established containing only pairs of species (e.g., Godoy et al. (2020)), we mostly measured polycultures with several different neighbour species (e.g., T. officinale surrounded by C. biennis,C. jacea and R. acetosa), as the experiment did not contain two species plots. However, we were able to estimate all pairwise interspecific interactions thanks to the modelling approach applied (see below). This approach is part of the so-called individual fitness models, which are commonly used to estimate interactions between pairs from multispecies data (Bimler et al., 2018; Garc ´ ıa-Callejas et al., 2020; Martyn et al., 2021). The effect of higher-order interactions (Mayfield and Stouffer, 2017) on coexistence was beyond the scope of this article and we therefore did not estimate them in our models but our interaction coefficients average across any higher order interactions that might occur. We visually estimated the cover of the focal species inside the rings as well as the cover of the other species. The cover of the neighbouring species was estimated individually for the eight target 26 1.2. Materials and methods Figure 1.2: Sampling method through the PaNDiv experiment. (a) Picture of the field with detail of the 40 cm diameter ring used to measure the cover of the focal (marked with a stick and labelled) and neighbour plants. (b) Sampling design, using a natural gradient of neighbour density (x-axis) in monocultures and polycultures to estimate the intra- and inter-specific interaction coefficients, respectively. The y-axis represents a gradient of neighbour diversity for the polycultures. In this example, we have used Anthoxanthum odoratum (Ao) as focal (bold letters), but each one of the 8 species is sampled both as focal and neighbour. Note that asterisks indicate that the focal species may also be present as neighbour in polycultures, contributing to estimating the intra-specific coefficients together with the information from monocultures. 27 1. Nitrogen, pathogens, and multispecies coexistence species we selected, and we estimated the abundance of all other species grouped into three functional groups (herbs, grasses, and legumes). We did not explore the likelihood of coexistence of these groups but we determined their competitive effects on our eight target species (see statistical models below). We measured plant cover at the beginning of the growing season and in June, at peak biomass, before the first cut. This therefore represents the majority of aboveground growth for each plant in that year. We measured plant cover for logistical reasons, but we converted the cover values to biomass for the modelling. To do so, we harvested (and oven-dried, at 60 ºC for two days) the aboveground biomass of eight focal plants per species and treatment, to estimate a relationship between cover and the biomass of each focal species. We obtained reasonably good fits for all and each of the individual species (R 2 = 0.6 ± 0.2, mean ± SD) (Fig. S1.1). Finally, we repeated the whole procedure in 2018 on newly-selected focal neighbourhoods, which meant we measured, across the two years, the cover and species composition of 4,096 neighbourhoods (2,048 each year). Once we obtained data from the experiment, we visually checked for outliers. We took the decision to remove the 2018 data from one species, Rumex acetosa, from further analyses because it declined strongly in biomass from February to June during this year (Fig. S1.2), most likely because conditions were very dry during 2018 and R. acetosa had already begun to senesce by June. Modelling this strong decline would indicate that this particular species is not able to grow at the experimental site (i.e., environmental filtering), which is not the case in general. We then modelled the interactions between all pairs of species. We used a modified Beverton-Holt function to estimate model parameters (Beverton and Holt, 1957). This modification was applied to model population size as biomass and not as numbers of individuals because the clonal growth of several of these perennial plant species makes it impossible to distinguish separate individuals in the field. However, the model retains the mathematical behaviour of prior work that has shown this function can be applied to both annual and perennial plant species (Lanuza et al., 2018; Cardinaux et al., 2018), and is also suited to apply a structural approach (see details in Appendix 4 in Saavedra et al. (2017)). Since we marked individuals already growing in the field (we did not plant out individuals), we needed to account for the initial biomass of each focal. Focal individuals with greater initial biomass might experience less competition and grow more, as bigger individuals are less sensitive to competition (Rees, 2013). Therefore, we modelled the biomass response ( N ) of the focal species i in the year t+ 1 as follows: Ni,t+1 Ni,t =λi(1 + θiNi,t) 1 + Pn j=1(αij +Ni,t)Nj,t (1.1) where λi is the per capita growth rate measured as biomass of i in the absence of interactions and αij represents the effect that each neighbour species has on i . The parameters θi and Ωi reflect the effect of the initial size of species i on its intrinsic growth rate ( λi ) and its interactions with neighbour species ( αijαij ), respectively. We fitted eqn. 1.1 to empirical observations using maximum likelihood methods 28 1.2. Materials and methods ( optimx function within the “optimx” package (v. 2022-4.30, (Nash and Varadhan, 2011); method = “Nelder-Mead”). We estimated λi (bounded to be positive as negative intrinsic growth rates for species growing alone do not make sense given that all species can persist at the site) and αij for each treatment separately (control, N addition, fungicide, and combined). αij is not bounded and can be negative (competition) or positive (facilitation) because the structural approach can compute niche differences when both types of interactions are combined. We fitted three different models for each species to deal with initial biomass effects. We considered that intial biomass ( Ω ) had: 1) no effect on the other model parameters (i.e., the effect of size is negligible), 2) a common size effect across all species, or 3) a species-specific size effect. To identify the best model for each species and treatment, we used AIC (Burnham and Anderson, 2002). All the details about the three different specifications of the function, and their AIC estimates for each species, are included in Table S1.2. Note that this modelling procedure was also done separately for the years 2017 and 2018. To better understand the confidence of our results, we calculated the Hessian matrix with the ghgen function within the “optimx” package and transformed its diagonal to obtain the standard error for all model parameter coefficients (95% confidence), assuming that the gradient is well approximated by a quadratic function. Mean model estimates and standard errors for the αij and λi parameters under each treatment are provided in Table S1.3—S1.6. 1.2.3 Predicting the effects of N enrichment on species coexistence with the structural approach With the field parametrisation of the dynamic model describing the population dynamics of interacting species, we applied the structural approach to explore how coexistence mechanisms varied across the experimental treatments. We followed previous definitions of Saavedra et al. (2017) to compute metrics of structural niche differences (SND), structural fitness differences (SFD), and feasibility (i.e., an equilibrium point which is assumed to be globally stable following previous approaches (Zhao and Luo, 2010; Saavedra et al., 2017); thus, here we consider a feasible combination to be one in which all species coexist). We did this for species pairs and for combinations of three to six species. The great advantage of the structural approach is that we can evaluate the probability of coexistence for any number of species at the same time if we have information on their intrinsic growth rates and the full matrix of pairwise interactions. For all species combinations, coexistence conditions are the same: a species combination can coexist when SND are greater than SFD 1.1. Importantly, we can compare the variation in SND and SFD across treatments, but we cannot compare them across combinations of different numbers of species, i.e., SND in a three species combination (triplet) cannot be compared with the SND in a four-species combination. This is because SND and SFD inherently vary as more species are included in the combination, as indirect interactions between species contribute to SND and these vary with the number of species included. We also evaluated the importance of indirect interactions (e.g., intransitivity and other types of species interaction chains) in 29 1. Nitrogen, pathogens, and multispecies coexistence allowing species to coexist using the metric called community-pair differential (Saavedra et al., 2017). Community-pair differential varies between − 1 and 1. A positive value indicates that indirect interactions promote coexistence, i.e., there are greater opportunities for coexistence in the multispecies combination than in the individual pairs that make up the combination. For example, in the case of an intransitive triplet, no individual pair could coexist, but the triplet could coexist. Negative values of the differential indicate the opposite, i.e., fewer opportunities for coexistence in the multispecies combination. These computations were performed for the mean model parameters, the lower and upper tails of the 95% confidence intervals, and for the years 2017 and 2018 separately. We finally used generalised linear models (GLMs) to evaluate the effect of N addition and fungicide application on these coexistence metrics, for each number of species in the combination. To determine what distribution to use for the residuals in the analysis of SND, SFD, community-pair differential, and feasibility metrics, we used the function descdist within the package “fitdistrplus” (Delignette-Muller and Dutang, 2015), which led us to fit a GLM with gamma errors (link = “identity”) for SND, SFD, and community-pair differential metrics, and binomial errors for the feasibility metric. A summary of the results from the different models for all metrics is available in Table S1.7—S1.26. All computational analyses were performed in R version 4.2.1 (R Core Team, 2022). 30 1.3. Results 1.3 Results 1.3.1 Determinants of feasibility across treatments Our experimental treatments resulted in communities with a wide range of structural niche differences (SND) and structural fitness differences (SFD). Nevertheless, the variation was higher for SFD than for SND regardless of the number of species considered. Accordingly, most communities were predicted to be unfeasible and not able to coexist (Fig. 21.3). We did not observe a clear threshold between feasible and unfeasible assemblages in multispecies communities mainly because indirect interactions promoted species persistence for some species combinations and hindered it for others. For instance, feasible communities were those in which we found 1) intra-specific strongly exceeded inter-specific interactions, 2) quasineutral-like dynamics in which all interactions and growth rates were equivalent but intraspecific interactions were slightly stronger than inter-specific ones, and 3) intransitive dynamics. Conversely, examples of predicted unfeasible communities were those in which inter-specific interactions were stronger than intra-specific interactions (i.e., priority effects), and 2) high variation in intrinsic growth rates among species destabilised the effect of niche differences. 1.3.2 Effects of nitrogen enrichment and fungal pathogen suppression on structural niche and fitness differences The effect of nitrogen (N) addition, pathogen suppression, and their interaction on changing SND compared to the control varied depending on the number of interacting species (Fig. 1.4a). Within species pairs, N addition and pathogen suppression increased SND when combined but did not differ from control conditions when applied alone. N addition slightly decreased SND in three-species combinations; however, the combination of N addition and pathogen suppression increased SND. For assemblages of four, five, and six species the pattern was more consistent: N addition, pathogen suppression, and their combination tend to increase SND compared to control conditions. Separating each year individually, we observed in 2017 that N addition, pathogen suppression, and their combination decreased SND. This result differs from both years combined and such differences are more pronounced with a higher number of interacting species (Fig. S1.3a). Analysing 2018 alone showed similar results to both years combined (Fig. S1.3b). We also explored how our experimental treatments modified SND at the lower and upper limits (95% CI), which rendered smaller and less variable SND values than the mean SND for both years (Fig. S1.4). On the other hand, SFD did not vary between treatments and control conditions for combinations of two species whereas for assemblages of three species only the combination of N addition and pathogen suppression slightly increased SFD. Similar to the case of SND, we observed parallel trends from four to six-species assemblages. Specifically, N addition, pathogen suppression, and their combination strongly decreased SFD (Fig. 1.4b). The separated analyses of SFD in 2017 and 2018 showed contrasting results. In 2017, the reduction of SFD compared to control conditions is repeated for assemblages of four, five, and six species (Fig. 31 1. Nitrogen, pathogens, and multispecies coexistence Figure 1.3: Diagram describing the relationship between the structural niche differences (SND) and the structural fitness differences (SFD) for all three-species combinations across treatments in our study, together with some examples of how these mechanisms relate to the underlying network of interaction coefficients. Data points in black and white represent unfeasible and feasible triplets of species, respectively. The triangles visually represent the structural approach (Saavedra et al., 2017) as thoroughly described in Fig. 1.1 and the main text. As a reminder, coexistence is achieved when the feasibility domain or SND (green triangles, defined by the structure of species interactions ( α )) are capable of accommodating the vector of intrinsic growth rates ( λ ; orange dot). The distance between λ and the centroid of the feasibility domain ( × ) visually approximates the SFD. Low values of SND usually generate unfeasible combinations (a), although even some small feasibility domains can have shapes that may accommodate intermediate values of SFD (b) due to the structure of the interaction network. Moreover, quasi-neutrality (c) is a particular example of a feasible combination where both SND and SFD tend to zero. Intransitive dynamics such as rock-paper-scissors (d) may also occur, in which the multispecies combination is feasible but all the pairwise interactions are unfeasible. Finally, even assemblages with high niche dimensionality can be unfeasible in combination with specific intrinsic growth rates (e), although the most likely scenario is that these intrinsic growth rates can be accommodated and a feasible assemblage produced (f). We show the values of the matrix of species interactions ( α matrix) and the vector of intrinsic growth rates (λvector) that underlies each representation. 32 1.3. Results S1.5a). However, we observed the opposite pattern in 2018; namely, all treatments increased SFD in assemblages of four-, five- and six-species combinations (Fig. S1.5b). Despite the contrasting results, analysing both years individually confirms the general pattern that variation of SFD differences between treatments was highly dependent on the number of interacting species. Finally, the exploration of how our experimental treatments modified the lower and upper limits (95% CI) of SFD showed smaller and more variable values but comparatively consistent with the mean values for both years, particularly in the upper limit (Fig. S1.6). Figure 1.4: Structural niche (a) and fitness (b) differences for all combinations of two to six species, in the four experimental treatments (control, nitrogen, fungicide, and combined). Each point represents a value of SND or SFD of a specific species combination under a particular treatment. Asterisks above boxes indicate significant differences (p < 0.05) between the corresponding treatment and the control (Table S1.7—S1.8). Note that the axis scales differ between panels because results can be compared across treatments within each species combination but not between species combinations (see text). Results represented correspond to mean estimates of interaction coefficients and intrinsic growth rates when data from years 2017 and 2018 were combined. For mean values of each year separately, the same representation and statistical analyses can be found in Fig. S1.3 and Table S1.9—S1.10 (SND) and S1.5 and Table S1.11—S1.12 (SFD). For the lower and upper limits (95% CI) of both years combined, the same representation and statistical analyses can be seen in Fig. S1.4 and Table S1.13—S1.14 (lower limit) and in Fig. S1.6 and Table S1.15—S1.16 (upper limit). 33 1. Nitrogen, pathogens, and multispecies coexistence 1.3.3 The role of indirect effects Both direct pairwise interactions and multispecies indirect interactions can affect the structural niche differences (SND), and the relative contribution of each interaction type can be evaluated by the community-pair differential. Overall, we observed small values of community-pair differential (i.e., on average across treatments, only 15% of the SND were due to indirect interactions, for three species combinations), and these community-pair differential values tended to zero with increasing species richness. Contrary to our initial expectation, this result suggests that indirect interactions have a minor role in promoting SND as more species are added (Fig. 1.5). Equally surprising is that we did not observe differences when specific treatments were compared. When decomposing the main results into 2017 and 2018, we observed the same overall decay in the importance of indirect effects with the number of interacting species and experimental treatments had no effect (Fig. S1.7). The exploration of the lower and upper limits (95% CI) of community-pair differential values corroborated the low importance of indirect interactions on promoting species coexistence because all values tended to zero in the lower limit regardless of the number of species and treatment considered (Fig. S1.8a), and upper limits reproduced the decaying importance of indirect interactions for promoting SND as the number of species increased (Fig. S1.8b). Figure 1.5: Community-pair differential for combinations from three to six species simultaneously, in the four experimental treatments (control, nitrogen, fungicide, and combined). Each point represents a value of community-pair differential of a specific species combination under a particular treatment. Asterisks above boxes indicate significant differences (p < 0.05) between the corresponding treatment and the control (Table S1.17). Results represented correspond to mean estimates of interaction coefficients and intrinsic growth rates when data from years 2017 and 2018 were combined. For mean values of each year separately, the same representation and statistical analyses can be found in Fig. S1.7 and Table S1.18—S1.19. For the lower and upper limits (95% CI) of both years combined, the same representation and statistical analyses can be seen in Fig. S1.8 and Table S1.20—S1.21. 34 1.3. Results 1.3.4 Consequences for diversity maintenance and changes in species composition The combination of SND and SFD resulted in a reduced proportion of feasible species combinations as the number of interacting species increased (Fig. 1.6). This result is particularly apparent for the control treatment, which does not sustain any feasible combination of six species. However, the nitrogen and fungicide treatments increased the number of feasible combinations when four or more species were considered. When we separated both years, we observed a more abrupt decay in the proportion of feasible combinations with the number of interacting species, particularly in 2018 when there are no feasible combinations of five or six species (Fig. S1.9). In 2017, the fungicide treatment seems to support the most feasible combinations, especially in combinations of three and four species. The exploration of the lower and upper limits (95% CI) in both years combined offered a substantial decrease in feasibility with almost no predictions of feasible combinations, especially for multispecies assemblages (Fig. S1.10). Figure 1.6: Feasible combinations (% of the total) of two to six species, in the four experimental treatments (control, nitrogen, fungicide, and combined). A feasible combination is one in which all species are predicted to coexist. Asterisks indicate significant differences (p < 0.05) between the treatments and the control (Table S1.22). The total number of possible combinations (combos) of two-six species is shown at the bottom. Results represented correspond to mean estimates of interaction coefficients and intrinsic growth rates when data from years 2017 and 2018 were combined. For mean values of each year separately, the same representation and statistical analyses can be found in Fig. S1.9 and Table S1.23—S1.24. For the lower and upper limits (95% CI) of both years combined, the same representation and statistical analyses can be seen in Fig. S1.10 and Table S1.25—S1.26. Together with the impact on the overall numbers of species predicted to coexist, our treatments did predict strong compositional changes. For instance, only five species were present in all species combinations under control conditions, with T. officinale and A. odoratum being unable to coexist in five- and six-species combinations and P. media in six-species combinations (Fig. 1.7a). Conversely, most species were evenly distributed across all feasible species combinations 35 Loss of species richness with land use intensity is explained by a reduction in niche differences 2 43 2. Multispecies coexistence in a gradient of land use intensity 44 Abstract Increases in land use intensity (LUI) reduce species richness. However, we have a poor understanding of how underlying coexistence mechanisms are altered by land use and whether diversity loss occurs due to changes in plant-plant interactions (competition and facilitation) or in species intrinsic growth rates. We expect that LUI could reduce stabilising niche differences and the indirect interactions that promote coexistence (e.g., intransitivity), while increasing competitive inequalities between species. To test the importance of these different processes, we use 8-yr time series from 150 grasslands differing in LUI to evaluate the role of direct and indirect interactions in promoting coexistence between 50 plant species. We show that LUI reduces the number of coexisting species mostly by causing a non-linear reduction in niche differences, rather than by enhancing competitive inequalities. However, surprisingly, niche differences remained important in stabilising coexistence between those species remaining at high LUI. Indirect interactions were generally less important than direct ones, and played a moderate role in promoting coexistence in smaller assemblages of species at intermediate LUI. Our models could accurately reproduce the decline in diversity seen with LUI, indicating that our time series approach captures the important interactions between species. By analyzing land use effects through recent advances in structural stability applied to community ecology we provide a more mechanistic understanding of its effects. Our results highlight the importance of identifying the niche differences that are lost with increasing LUI, to better predict and manage effects of land use on biodiversity. Keywords Global change, land use intensity, biodiversity loss, multispecies coexistence, structural stability 45 2. Multispecies coexistence in a gradient of land use intensity 46 2.1. Introduction 2.1 Introduction Changes in land use intensity (LUI), for instance increases in fertilisation, mowing (e.g., early and frequent cuts), or grazing intensity, are often a major threat to biodiversity (Rockstr ¨ om et al., 2009). Although many studies have shown that plant diversity declines (Newbold et al., 2015; Allan et al., 2014; Crawley et al., 2005) and species composition changes (Pakeman, 2011; Dormann et al., 2007) as land use is intensified, we have a poor understanding of the mechanisms responsible, which hampers efforts to predict, conserve, and manage land use effects. Several studies have tried to understand mechanisms of community assembly from changes in functional diversity and have shown that functional diversity declines with land use intensification as communities at high land use intensity become dominated by functionally similar species (Pakeman, 2011; Mayfield et al., 2010; Lalibert ´ e and Tylianakis, 2012). For instance, the reduction of plant canopy height, leaf dry matter content, and seed mass associated with land use intensification has been interpreted as an environmental filtering process (Pakeman, 2011). However, it can be challenging to infer precise coexistence mechanisms from functional trait distributions (Mayfield and Levine, 2010). An alternative perspective is to take a demographic approach in which changes in population abundance over time are analysed to estimate the role of species interactions and environmental conditions in determining species coexistence, i.e., the probability that interacting species can persist over time (Hallett et al., 2019; Song et al., 2020b; Wainwright et al., 2018). In general, land use intensification could alter species intrinsic growth rates, and the strength of competitive and facilitative interactions (Olsen et al., 2016; Adler et al., 2018; HilleRisLambers et al., 2012) between them, by changing the biotic and abiotic environment. If the altered environmental conditions in intensively managed grasslands result in negative intrinsic growth rates for some species, theory predicts that they will become extinct, with rare species more likely to be lost by chance (Pakeman, 2011; Hallett et al., 2019; Selwood et al., 2015). However, land use intensification can cause further diversity loss by altering species interactions and the balance of intrinsic growth rates between them. Land use intensification can therefore disrupt two main mechanisms of coexistence (Saavedra et al., 2017): 1) stabilising mechanisms, defined as niche differences, stabilise the population dynamics of interacting species, by causing negative frequency-dependent growth; and 2) equalising mechanisms, defined as fitness differences, drive competitive similarity among species when differences in intrinsic growth rates are reduced. Various aspects of land use intensification may alter niche and fitness differences. Some well documented cases are when fertilisation favors fast-growing species, and thereby increases competitive inequalities, by increasing the intrinsic growth of good light competitors at the expense of other species (i.e., LUI could increase fitness differences) (DeMalach et al., 2017). These increases in competitive asymmetry make it harder for species to coexist, unless superior competitors suffer self-limitation, which occurs when intraspecific exceeds interspecific competition (i.e., niche differences stabilise competition between species). All processes that cause intraspecific competition to exceed interspecific (e.g., specialised resource uptake, natural enemies, etc.) contribute to these niche differences. Such 47 2. Multispecies coexistence in a gradient of land use intensity self-limiting niche differences could also be reduced in intensively-managed grasslands: for instance, fertilising grasslands could prevent trade-offs in nutrient use from stabilising coexistence, by reducing the number of resources, such as N, P, and K, that plants compete for (Wilson and Tilman, 1993). These changes would lead to smaller niche differences between plant species, and therefore, less stable coexistence between them (Adler et al., 2018; Saavedra et al., 2017). On the other hand, diversity loss can lead to an increase in specialist pests and pathogens (Bagchi et al., 2014), which might stabilise coexistence between remaining species. Many other changes in coexistence mechanisms are possible and land use intensification could therefore lead to complex changes in multiple mechanisms stabilising the dynamics of interacting species. However, this range of different mechanisms can be collectively explored by quantifying overall niche differences along gradients of land use intensity. In multispecies communities, species can also coexist through two types of indirect interactions. The first possibility is intransitive competition, e.g., the rock-paper-scissors game in which species A excludes B, species B excludes C, and species C excludes A, resulting in no universally superior competitor (Soliveres et al., 2015). Weaker forms of intransitivity are also possible and could promote coexistence in combination with niche differences (Gallien et al., 2017; Godoy et al., 2017). The second possibility is that a superior competitor differentially harms two inferior competitors and thereby indirectly facilitates one of them (Miller, 1994). For instance, species B and C cannot coexist together alone but if species A harms species B more than C, then it can indirectly facilitate species C and allow the triplet to coexist (Miller, 1994). Land use may reduce the importance of these indirect interactions for coexistence (Soliveres et al., 2015; Laird and Schamp, 2006): for example, land use intensification can lead to more hierarchical competitive networks (DeMalach et al., 2017; Wilson and Tilman, 1993), in which a few species dominate and exclude the others, because it leads to greater differences in individual size and therefore a greater imbalance in competitive effects (Rees, 2013). In addition, land use could disrupt other types of indirect interactions because it reduces facilitative interactions between species (Mayfield et al., 2010; Olsen et al., 2016). In sum, land use intensification is expected to reduce coexistence via three main mechanisms: 1) enhancing competitive inequalities between species, 2) destabilising coexistence by reducing niche differences, and/or 3) reducing the opportunity for indirect interactions to promote species coexistence. No studies have been able to determine the relative importance of these different mechanisms in driving the species loss (i.e., reduction in number of species) frequently observed with land use intensification. A key objective is therefore to better understand how coexistence mechanisms emerging at the pairwise and at the multispecies level vary simultaneously along a LUI gradient and how this contributes to species loss in managed grasslands. Here, we address this gap by applying a structural stability approach (Saavedra et al., 2017; Rohr et al., 2014) to a unique dataset from the German Biodiversity Exploratories. This includes data on changes in abundance/cover for 50 perennial plant species, across 8 years, in 150 grasslands differing in their degree of land-use intensification (Bl ¨ uthgen et al., 2012; Fischer et al., 2010). In these grasslands, 48 2.1. Introduction we observe strong declines in species richness with increases in LUI (Newbold et al., 2015; Flynn et al., 2009) (Fig. S2.1), similar to the decline in richness seen in a long term fertilisation (Crawley et al., 2005). The structural approach allows us to estimate the three types of mechanism that contribute to coexistence, for assemblages of any number of species (see Fig. 0.1 for a detailed description). Briefly, both niche differences and indirect interactions jointly determine the size of the coexistence region (i.e., the feasible domain), and the larger the region, the more likely species are to coexist. At the same time larger differences in the ability of a species to increase its population (i.e., differences in intrinsic growth rate) make coexistence more difficult (Saavedra et al., 2017; Godoy et al., 2018). We first use the changes in species abundances (i.e., proportion of cover) over time to estimate species intrinsic growth rates, and all pairwise interactions between species, and to assess how LUI changes these growth rates and pairwise interactions. We then couple a population model with the structural stability approach to estimate how land use changes coexistence mechanisms. We do this by calculating the expected structural niche differences, structural fitness differences and indirect interactions for all sets of two and three species, at multiple points along the LUI gradient. To validate our models and to show the importance of these modelling tools for practitioners and conservation biologists, we use them to test whether the predicted number of species that is expected to coexist along the land use intensity gradient is similar to the observed pattern in the decline of number of species with LUI (Newbold et al., 2015; Flynn et al., 2009) (Fig. S2.1). 49 2. Multispecies coexistence in a gradient of land use intensity 2.2 Materials and methods 2.2.1 Study system The Biodiversity Exploratories (www.biodiversity-exploratories.de) project (Fischer et al., 2010) is a research initiative that has established 150 permanent grasslands plots of 50 × 50 m in three different regions across Germany: the UNESCO Biosphere Area Schw ¨ abische Alb (south-west), Hainich National Park (central) and the UNESCO Biosphere Reserve Schorfheide-Chorin (north-east). All three regions have a similar climate (a range of 3 ºC in mean annual temperature, and 500 to 1000 mm annual precipitation—more details provided in Fischer et al. (2010)). These 150 plots have been managed as grasslands for at least 20 years prior to the start of the project. Farmers and landowners provide information on the intensity of land management activities, including fertilisation, mowing, and grazing. These sites, which range from seminatural to intensively managed grasslands, are either mown or grazed. Grazed plots have cattle, horses, or sheep, with varying numbers of animals and duration of grazing. Grazing intensity was quantified as the number of livestock units ha −1 year −1 . Mown plots are cut one to three times per year. Some grasslands are also fertilised and the intensity is quantified as the amount of organic and inorganic nitrogen added to the grassland (Bl¨ uthgen et al., 2012). With this detailed land use information, prior work has quantified a compound index of land use intensity (LUI), which integrates the intensity of fertilisation (F), the mowing frequency (M), and the intensity of grazing (G) for each grassland plot (Bl ¨ uthgen et al., 2012; Fischer et al., 2010). For each plot, an individual LUI component (F, M, or G) was standardised relative to its mean across all three regions and all years considered. The compound LUI is the sum of the three standardised components (more details are provided in Allan et al. (2014)). The main advantage of this index is to summarise different, correlated aspects of land use into a single metric. We use this index to estimate the effect of human mediated actions on species interaction coefficients and intrinsic growth rates (see next section, estimation of species interactions). The minimum LUI of 0.5 could be produced by mowing every 2 y, fertilising at the rate of 6 kg of N ha −1 y −1 , or grazing one cow ( > 2 y old) per hectare for 30 d (or one sheep per hectare for the whole year). An intermediate LUI of 1.5 would equate to around two cuts per year, the addition of 60 kg of N ha −1 y −1 , or grazing one cow per hectare for most of the year (300 d). A high LUI of 3.0 could be produced by grazing by three cows per hectare for most of the year (300 d) and fertilising at the rate of 50 kg of N ha −1 y −1 or by cutting three times and fertilising with 130 kg of N ha −1 y −1 . For more details see Allan et al. (2014); Fischer et al. (2010). Within the period of 8 years analysed, LUI varied in our study between 0.5 to 3.0 across plots (Fig. S2.1). 2.2.2 Estimation of species interactions The main aim of our study was to investigate how coexistence mechanisms, at local scales, vary along the LUI gradient. To achieve that aim, we need to estimate species intrinsic growth rates, the matrix of species interactions, including intra- 50 2.2. Materials and methods and inter-specific interactions, and the effects of LUI on both species interactions and intrinsic growth rates. By doing a space for time substitution, we analysed a time-series of changes in the proportional cover of the 50 most common species (i.e., the most frequently observed species across all plots) between 2008 and 2015. Percentage cover (converted to proportions) was measured for each species and plot in 4m × 4m subplots at peak biomass by (Socher et al., 2012). To obtain estimates of 1) species intrinsic growth rates, 2) pairwise species interaction coefficients, and 3) the effect of LUI on these parameters, we fitted generalised linear mixed models (GLMM). Using a similar approach to (Tredennick et al., 2017; Clark et al., 2020; Grenfell et al., 1998), we analysed the log ratio of cover at time t + 1 and cover at time t for each focal species (as the dependent variable). The log ratio of cover was modelled as a linear function of LUI, together with the cover of each species in the community (including itself) at time t , and we further included the interaction between LUI and cover of all species. The model also included plot nested within region as a random effect and accounted for temporal autocorrelation using an ARMA structure corAR 1 ( t + 1). We also tested for autoregressive components with greater time lags ( t + 2 and t + 3) but they resulted in a worse statistical fit. Models were fitted using the function lme in the R package ”nlme” (Pinheiro, 2011). The general form of the statistical model fitted is as follows: ln Coveri,t+1 Coveri,t ∼ri,LUI + n X j=1 αij,LUI ×Coverj,t,random =plot/region+corAR1(year) (2.1) In this statistical model, the intercept ( ri,LUI ) was interpreted as the intrinsic growth rate because it estimates the maximum ability of a species to change in cover between two years ( t +1 /t ), at different LUIs, in the absence of any neighbours in the previous year ( t ), i.e., at intra- and inter-specific cover of 0. Because we fitted a separate model for each species i , values of the intercept can vary across species and LUI. In turn, the alpha coefficients ( αij,LUI ) describe how positive (facilitation) or negative (competition) per-capita interactions with conspecific and heterospecific neighbours change with plant cover over time. Again, because our model is fitted for each species i independently, it includes the possibility that LUI differently modifies each pairwise interaction coefficient and therefore changes the whole network of species interactions at the community level. Values of species’ intrinsic growth rates at different LUI values ( ri,LUI ) are inferred from the regression and lack of data did not allow us to corroborate them empirically. However, we did test the sensitivity of our results to uncertainty in the parameter estimation by calculating coexistence mechanisms (i.e., structural niche and fitness differences) using parameter values at the 2.5% and 97.5% confidence intervals around each parameter (Figs. S2.5 and S2.6 and Table S2.2 and S2.3). We observed that LUI linearly changed intrinsic growth rates for all species. In fact, some intrinsic growth rates were negative at certain LUI values indicating that LUI can act as an environmental filter excluding certain species (Kraft et al., 2015a). Finally, we did not explicitly model dispersal in our study because there is no information available to assess the influence of propagule pressure on species 51 2. Multispecies coexistence in a gradient of land use intensity novel opportunities for species coexistence, we computed a pair of related metrics for sets of three species (Saavedra et al., 2017; Song et al., 2018). These metrics are community-pair differential and community-pair overlap. Briefly, communitypair differential is a metric ranging between -1 and 1 that quantifies whether the coexistence region increases or decreases with indirect interactions, i.e., whether indirect interactions promote or reduce coexistence. Community-pair overlap is a metric that quantifies the proportion of the coexistence region for a triplet that is similar to its constituent pairs. It ranges between 0 when triplets and pairs are completely different in their coexistence region, i.e., they coexist in different ways, and 1 when regions are exactly the same. Neither of the metrics quantifying the effect of indirect interactions on coexistence varied along the LUI gradient in sets of three species, indicating no change in the importance of indirect interactions. Community-pair differential values (which quantify whether indirect interactions decrease or increase the coexistence region) were close to zero (Fig. 2.2a), indicating that indirect interactions do not enhance coexistence. In addition, community-pair overlap values (the proportion of the coexistence region for a triplet that is similar to its constituent pairs) ranged from 0.5 to 0.8 (Fig. 2.2b), indicating that the majority of three species combinations can coexist because each of the three pairs that compose the triplet can also do so. These results contrast with prior work done in this system, which suggested a major role for indirect interactions, and particularly intransitivity (Soliveres et al., 2015), in promoting species richness. One potential explanation for these contrasting results is that different indirect interactions (e.g., intransitivity against indirect facilitation) have contrasting effects on coexistence and cancel out to reveal no overall effect of indirect interactions in this analysis. Another possibility could be that the current analysis included subdominant species, while the previous analysis (Soliveres et al., 2015) focused only on intransitivity between dominant species. As intransitive competition has been shown to be stronger between dominant species (Soliveres et al., 2015) this might explain the low importance of indirect interactions in the current study. Further work is needed to explore a wider range of indirect interactions and to determine how they change with land use intensification. It is also possible that indirect interactions become more important for coexistence as we consider larger sets of species. This is because the likelihood of finding a configuration of species interactions that maintains biodiversity is higher in larger sets of species. We therefore also computed community-pair differential and community-pair overlap metrics for assemblages of 5, 7, and 11 species, but only focused on the 26 most abundant species in this analysis to save computing power. We selected these numbers because it has been shown that indirect interactions promote coexistence when there is an odd number of species in the community, under the assumption of a random tournament (Gallien et al., 2017; Allesina and Levine, 2011), and 11 species is on average the maximum number of species found per plot (from the set of 26 species analysed). We found that community-pair differential remained low in combinations of 5, 7, and 11 species. This suggests that indirect interactions did not increase or decrease the opportunities for coexistence amongst larger sets of species (Fig. 2.2a). However, indirect interactions did change the shape of the feasibility domain, particularly 58 2.3. Results and Discussion Figure 2.2: Community-pair differential (i.e., the importance of indirect effects for coexistence, panel a) and community pair overlap (i.e., the importance of pairwise only interactions against indirect effects for coexistence, panel b) across the LUI gradient for assemblages of 3, 5, 7, and 11 species. The colors show species assemblages of different sizes. Values correspond to the mean for each number of species. Standard errors are not provided to facilitate visualisation. 59 2. Multispecies coexistence in a gradient of land use intensity for larger assemblages at intermediate LUI (community-pair overlap at LUI = 1.5 was ∼ 0.20 and 0.05 for 7 and 11 species respectively, i.e., 80% and 95% of the coexistence region due to indirect interactions do not overlap with the coexistence region promoted by direct interactions) (Fig. 2.2b). This indicates that indirect interactions that occur in multispecies assemblages change the way in which species coexist but they do not lead to more or less coexistence, perhaps because different indirect interactions have opposing effects on coexistence. Overall, these results show that direct pairwise interactions are the main mechanism promoting coexistence across assemblages; however, the opportunities for coexistence that indirect interactions create are completely different compared to those created by direct interactions, particularly at intermediate LUI. We then analysed whether we could reproduce the observed decline in species richness with LUI, based on our estimates of how species interactions varied. This allows us to validate our estimates of species interactions by using them to predict diversity. As we found that pairwise interactions were most important for coexistence overall, we focused on the coexistence predicted between species pairs (see the slope of the decline of species richness with LUI; Fig. S2.1). We computed the cumulative number of species predicted to coexist in any pairwise combination. Figure 2.3: The average number of species (grey line) ( ± standard deviation) observed to cooccur in a plot, the maximum number of species predicted to coexist according to the structural stability approach (in blue), and the maximum number of species predicted to coexist according to null expectations by randomising species interactions (in orange). Predictions from the models parameterising species interactions fit well with the observed decline in species richness with LUI. However, randomising interactions between species leads to much higher coexistence than observed. See methods for details on how the predictions were made. 60 2.3. Results and Discussion For instance, if at LUI = 0.5, species A coexists with B, B excludes C, but species C coexists separately with D, we assume that all four species can coexist within the grassland. We then repeated this procedure along the LUI gradient (from 0.5 to 3.0). This approach has been recently shown to reproduce the increase of species richness with area in spatially variable environments because it implicitly allows different species pairs to coexist in different patches (Garc ´ ıa-Callejas et al., 2021). Despite the simplicity of the model, we were able to reproduce the decline in species richness with LUI (Fig. 2.3). Note that we only calculated species richness based on the 50 most common species analysed in our study, which resulted in an initial slight increase in observed species richness with LUI, followed by a strong decline. However, including all species results in an exponential decline in species richness with LUI because many rare species are only present at the lowest LUI (Allan et al., 2014). It was not possible to include these rare species in our analysis as we lacked data to estimate changes in their abundance. Further experimental work would be extremely valuable in identifying mechanisms of coexistence between rare species. Overall, we found that the cumulative number of species predicted to coexist across the LUI gradient was close to the observed species richness, and different from randomisation predictions which did not capture the reduction in species richness with LUI. Instead the random expectation would be for variation in diversity but no net change with LUI. The only deviation was that we predicted slightly fewer coexisting species than observed at intermediate to high LUI values (LUI = 2.25—2.50) (Fig. 2.3). This difference is most likely due to the reduction in intraspecific competition at intermediate LUI (Fig. S2.2). Although the comparison between predicted and observed diversity suggests that we might have slightly underestimated the strength of niche differences at intermediate to high LUI, our overall ability to describe the decline in species richness with LUI provides strong evidence that the theoretical approach used here can well describe the observed patterns. By using time series data from 150 plant communities, we were able to estimate how land use intensification causes declines in plant diversity. Our findings provide strong evidence that a decline in niche differences with LUI explains the loss of diversity and, surprisingly, that changes in competitive ability play a relatively minor role. Determining what niche differences are lost as we intensify grassland management is a crucial next step that could enable us to offset the negative effects of land use on diversity. We also find that niche differences and interspecific facilitation maintain coexistence even in heavily modified ecosystems. Very little work has been done on how species coexist at high land use intensity, but our results show that, although diversity is reduced, coexistence between the remaining species is highly stable. Uncovering the processes that drive niche differences at high LUI might allow us to better understand and even optimise the ecosystem functioning of intensively managed grasslands. Our coexistence models were able to reproduce the observed diversity declines suggesting our results are robust. The approach of using time series to estimate species interactions and coexistence mechanisms could therefore be widely applied to derive a more mechanistic understanding of global change effects on biodiversity. 61 2. Multispecies coexistence in a gradient of land use intensity Associated publication • Godoy O, Granjel RR, van der Plas F, Soliveres S, Penone C, Saiz H, H ¨ olzel N, Prati D, Fischer M, Allan E. 2022. Loss of species richness with land use intensity is explained by a reduction in niche differences. bioRxiv: 2022.12.15.520291. Open science • R code: www.github.com/oscargodoy/LUI coexistence • Plant cover data: www.bexis.uni-jena.de/ddm/data/Showdata/27386 • LUI data: Vogt et al. (2019) Author contributions O.G., F.v.P., S.S., and E.A. conceptualised and formalised the idea. O.G., R.R.G., and E.A. designed the methodology and carried out the investigation. O.G. and R.R.G. computed the analyses and prepared the visualisation of results. O.G. and E.A. wrote the original draft, with subsequent review and editing from R.R.G. and all authors. Acknowledgements Special thanks to Serguei Saavedra for his guidance with the correspondence between the population and the statistical models, to the “Centro Inform ´ atico Cient ´ ıfico de Andaluc ´ ıa” (CICA) for the High-Performance Computing service and to Marcial Escudero for his assistance in using high-performance computation. Gratitude to all members of the Biodiversity Exploratories. Additional funding O.G. acknowledges financial support provided by the Spanish Ministry of Economy and Competitiveness (MINECO) and by the European Social Fund through the Ram ´ on y Cajal Program (RYC-2017-23666). S.S. was supported by the Spanish Government under a Ram ´ on y Cajal contract (RYC-2016-20604). H.S. is supported by a Mar ´ ıa Zambrano fellowship funded by the Spanish Ministry of Universities and European Union-Next Generation plan. This work has been partly funded by the DFG Priority Program 1374 “Biodiversity-Exploratories” (DFG-Refno. Po362/18-3). Supplementary material Supplementary figures and tables related to this chapter are available in Appendix 2. 62 Self-limitation and facilitation drive plant coexistence in diverse communities, outweighing herbivory and higher-order interactions 3 63 3. Self-limitation and facilitation in a complex herbivore-plant network 64 Abstract Highly-diverse ecological communities, such as plant-herbivore systems, involve a wide range of interactions that vary in sign, strength, and species. Due to the complexity of describing these interactions, it is difficult to understand which specific pairwise (plant-plant and plant-grasshopper) and higher-order interactions are critical in determining the long-term viability of plants (feasibility). In this study, we applied a regularisation approach to a detailed field experiment to identify the strength and sign of specific interactions among 36 plant species and 6 grasshoppers, which describe the spatiotemporal variation in plant cover. We then used this information along with theoretically informed metrics (structural niche and fitness differences) to understand what type and structure of species interactions drive the determinants of community feasibility. We found that only a subset of interactions across all types was selected, and among them, plantplant interactions played the predominant role in promoting the feasibility of modules of three and four plant species in the community. Within this key plant interactions compartment, the opportunities for species to coexist increased with facilitation and stronger intraspecific effects. Our results suggest that the specific structure of plant-plant interactions is key in promoting the feasibility of diverse systems with plant-herbivore interactions. Keywords Group lasso regularisation, grasshoppers, interaction metrics, multispecies coexistence, perennial plants, species modules, structural stability 65 3. Self-limitation and facilitation in a complex herbivore-plant network 66 3.1. Introduction 3.1 Introduction Understanding the role of ecological interactions in promoting species persistence (i.e., long-term positive abundance) is a longstanding research topic in ecology (Chesson, 2000; Saavedra et al., 2017). Theory predicts that coexistence between two competing species can be achieved when they limit their peers more than their competitors (i.e., intraspecific competition is greater than intraspecific competition) (Chesson, 2000; Barab ´ as et al., 2016). However, natural ecological systems are not composed of two but of multiple species belonging to different trophic levels, forming complex interaction structures that challenge our ability to interpret the mechanisms of coexistence (May, 1972; Mougi and Kondoh, 2012; Levine et al., 2017). Ecologists have addressed this challenge through three main perspectives. The first one aims to analyse the prevalence of positive, neutral, or negative interactions in ecological communities and the effect of the strength of such interactions on species’ growth rates (e.g., McCann et al. (1998); Neutel et al. (2002); Losapio et al. (2021)). The second perspective aims to compare whether interactions that occur horizontally within species of the same trophic level (e.g., direct and indirect competition) are more relevant to maintain species diversity than vertical interactions that occur between species of different trophic levels (e.g., herbivory) (e.g., Godoy et al. (2018); Song et al. (2018); Bartomeus et al. (2021)). Finally, the last perspective examines the configuration of species interactions within modules across whole networks (i.e., a combination of n -species that interact within a more complex, m -species network, with n being any natural number from 2 to m ) (e.g., Meli ´ an et al. (2009); Stouffer and Bascompte (2010)). Although these three approaches have been critical in providing us with important insights into the drivers of species persistence, there are still several main challenges that need careful consideration. It is a common practice to define the sign of interactions aprioristically, based on preconceived notions, which can lead to incorrect assumptions about the nature of such interactions (Bastolla et al., 2009; Rohr et al., 2014). For example, it is commonly assumed that pollinators and herbivores exert either positive or negative effects on plant fitness, but this is not always the case. Pollinators can also hinder plant fitness if their densities are too high (Magrach et al., 2017), whereas herbivores can increase, rather than decrease, per capita growth rates of perennial plants by pruning leaves and stems, potentially boosting biomass production (Karban and Strauss, 1993; Meyer, 1998; Kim et al., 2013). Therefore, a more accurate and nuanced understanding of the role of the interaction sign in species persistence can be obtained by empirical estimations rather than by theoretical assumptions (Godoy et al., 2018). In this regard, a recent study by Losapio et al. (2021) showed that alpine plant communities harbour more biodiversity when common modules including both positive and negative interactions occur. Nonetheless, the distribution of interaction strengths within each module together with their sign is key to predicting species’ persistence (May, 1972; Barab ´ as, 2021). Previous knowledge, mostly from food webs, has shown that the persistence of diverse systems arises when many weak and only a few strong interactions occur (McCann et al., 1998; Berlow, 1999; Neutel et al., 2002). Complementary, diverse systems are also found to be stable if their interactions show low variability irre- 67 3. Self-limitation and facilitation in a complex herbivore-plant network scenario accounts only for plant interactions but included both pairwise and those selected higher-order interactions between plant species, the βα layer. For the third scenario, we focus only on pairwise interactions but combine plant-plant interactions together with grasshopper-plant interactions, the γ layer. In the fourth scenario, we test a combination of multitrophic pairwise effects (plant-plant and plant-grasshopper interactions) together with higher-order interactions among plants. Finally, in the fifth scenario, we investigate the full empirical community by incorporating into the fourth scenario the higher-order interactions caused by grasshoppers, the βγ layer. Note that these scenarios range from simplicity to complexity hierarchically (i.e., βγ cannot take place without γ interactions, but vice versa). 3.2.6 Metrics for understanding the configuration of species interactions within modules In the last step of the analyses, we explored how the configuration of species interactions drives SND and SFD, the structural mechanisms of plant coexistence. To do so, we computed four different metrics for each three-species interaction module. To describe the shape of the distribution of interactions within a module, we first calculated its skewness and kurtosis. On the one hand, skewness is a measure of the symmetry of a distribution. It describes the extent to which a distribution is asymmetrical or skewed to one side or the other. Positive skewness indicates that the distribution has a higher concentration of values on the left side, with fewer values on the right side, whereas negative skewness indicates the opposite. On the other hand, kurtosis is a measure of the peakedness or flatness of a distribution. It describes the extent to which a distribution is more (positive or excess kurtosis) or less (negative or platykurtic kurtosis) peaked than a normal distribution, meaning that values in the distribution are respectively more or less concentrated near the mean. Together, they can be used to interpret the prevalence of weak and strong effects in distributions of interaction coefficients. Provided that theory predicts that coexistence requires interaction modules in which intraspecific effects are greater than interspecific effects, we measured in addition as the third metric the strength of intraversus interspecific interactions using an index called diagonal dominance (dom), which can be formulated as follows: dom = log10 Pn i=1 |Ai,i| n Pn i=1 Pn j=1 |Ai,j|−Pn i=1 |Ai,i| n(n−1) ,(3.2) where A is any given square matrix of coefficients describing a species interaction module and n is the number of columns in such a matrix. This index takes the mean of all diagonal coefficients, in absolute value, and divides it by the mean of non-diagonal coefficients, in absolute value. Therefore, positive values of diagonal dominance indicate stronger intraspecific than interspecific interactions and vice versa, with most of the values comprehended between -3 and 3 in our system. 74 3.2. Materials and methods Finally, the fourth and last metric computed was the positive-negative ratio (PNR), which describes the relative strength of positive versus negative effects within an interaction matrix and can be formulated as follows: PNR = log10  n P i=1 [Ai|Ai>0] m P j=1hAj|Aj<0i  .(3.3) Here, negative values of PNR indicate greater negative effects and positive values indicate greater positive effects in the matrix. In our system, most values of PNR range from -3 to 3. Thus, this metric can be used to assess the contribution of mutualisms and antagonisms across multispecies interaction modules in our empirical network. These four metrics were computed for all modules of three and four species under the full species interaction network (scenario 5). We then assessed with generalised linear models (GLMs) their individual and joint contributions to promoting SND and SFD. To determine the distribution of the residuals we used the function descdist within the package “fitdistrplus” (Delignette-Muller and Dutang, 2015), leading us to fit a GLM with gamma errors (link = “log”) for both SND and SFD. A summary of the results of these analyses is available in Table S3.2-S3.3, and specific results for modules of four species are available in Table S3.4-S3.5. All computational analyses were performed in R version 4.2.1 (R Core Team, 2022). 75 3. Self-limitation and facilitation in a complex herbivore-plant network 3.3 Results Figure 3.2: Components of our empirical system and their associated coexistence outcomes. (a) Representation of the inferred empirical interaction network for the full community. Cyan circles represent the set of 35 plant species and coral circles stand for the six grasshopper species. Notice that interactions between species can be either positive (solid links) or negative (dashed links) and that the width of such links is directly proportional to the strength of the effect. (b) Diagram showing all the different interaction types occurring within our empirical network: a plant-plant interaction α layer composed of intra-specific ( αintra ) and inter-specific ( αinter ) interactions, a γ layer accounting for the effects of grasshoppers on plants, and a β layer corresponding to the higher-order effects on plant-plant interactions. (c) Histograms show the density of positive and negative interactions for each of the aforementioned interaction layers. As zero values have been removed from the visualisation, we show a connectance value (Cn) for each layer that indicates the proportion of values different from zero. (d) Mechanisms of coexistence (structural niche and fitness differences) and predicted coexistence outcomes (feasibility) for the full community, i.e., including all interaction layers, calculated for all plant triplets (i.e., combinations of three plant species). 76 3.3. Results The regularisation approach allowed us to identify the structure of biotic interactions that characterise the complex network of relationships between the plants and grasshoppers in our experimental system (Fig. 3.2a). This characterisation connects the selected biotic interactions with the growth and population dynamics of plant species. Across types of interactions, including plant-plant pairwise interactions or α , pairwise grasshoppers effects on plants or γ , and the aggregated higher-order interactions (HOIs) of plants and grasshoppers or β ; we observe similar levels of connectance (i.e., percentage of interactions with a net effect different from zero). For instance, the connectance value of the α layer is Cn = 0.4, which means that 60% of interactions had a net effect of zero) (Fig. 3.2b). However, it is worth noting that the connectance showed for β is achieved by the joint effect of 36 different interactions matrices (30 matrices of plant HOIs and 6 matrices of grasshopper HOIs), which individually show very low connectance (Cn = 0.043 ±0.014). Figure 3.3: Contribution to the coexistence mechanisms ((a) structural niche differences and (b) structural fitness differences) of the effects of grasshoppers on plants ( γ ) and the additive higher-order effects of plants on plant-plant interactions ( β ). This contribution is measured as the change (%) in the coexistence mechanisms compared to a community with only plant-plant interactions (α). Whiskers represent the standard error for each bar. Regardless of the interaction type, all of them displayed both positive and negative signs but their distribution of strengths varied. While we observed 77 3. Self-limitation and facilitation in a complex herbivore-plant network a wider spread for plant-plant ( α ) (range from -0.57 to 0.59) and grasshopperplant interactions ( γ ) (range from -0.32 to 0.47), the distribution of HOIs ( β ) was concentrated around zero, meaning that the majority of effects were very weak. With this parameterisation, we explored how the different types of interactions drive the structural niche and fitness differences (SND and SFD) and determine the feasibility of three plant species modules (i.e., triplets). Following theoretical expectations, we observed that the majority of triplets predicted to be feasible followed the condition that SND overcomes SFD (Fig. 3.2c). In addition some triplets were predicted to be feasible despite showing low SND and high SFD. This was likely due to the positive effect of indirect interactions such as intransitivity (rock-paper-scissors dynamics) on promoting feasible triplets. Finally, it is worth Figure 3.4: Structural niche differences (a—f) as a function of four metrics (skewness, kurtosis, diagonal dominance, and PNR) measuring the configuration of interaction strengths within each three-species module in the community. Skewness describes how much a distribution is shifted to the left (stronger negative interactions) or right (stronger positive interactions) of its centre. Kurtosis describes the peak and tails of a distribution: higher values indicate a higher peak and thinner tails (more weak interactions) and lower values indicate a lower peak and thicker tails (more strong interactions). Diagonal dominance is an index describing the relative difference in the strength of the diagonal of the interaction matrix (intraspecific interactions) compared to its non-diagonal elements (interspecific interactions), with logarithmic values below zero indicating higher intrathan interspecific effects and vice versa. Positive-negative ratio (PNR) is an index describing the preeminence of positive (PNR > 0) or negative (PNR < 0) interactions within each module. See methods for a description of each metric. 78 3.3. Results noting that we observed many triplets not predicted to be feasible and these correspond to a range of situations in which low SND do not overcome low SFD as well as high SND do not overcome SFD (Fig. 3.2c). Figure 3.5: Structural fitness differences (a—f) as a function of four metrics (skewness, kurtosis, diagonal dominance, and PNR) measuring the configuration of interaction strengths within each three-species module in the community. Skewness describes how much a distribution is shifted to the left (stronger negative interactions) or right (stronger positive interactions) of its centre. Kurtosis describes the peak and tails of a distribution: higher values indicate a higher peak and thinner tails (more weak interactions) and lower values indicate a lower peak and thicker tails (more strong interactions). Diagonal dominance is an index describing the relative difference in the strength of the diagonal of the interaction matrix (intraspecific interactions) compared to its non-diagonal elements (interspecific interactions), with logarithmic values below zero indicating higher intrathan interspecific effects and vice versa. Positive-negative ratio (PNR) is an index describing the preeminence of positive (PNR > 0) or negative (PNR < 0) interactions within each module. See methods for a description of each metric. We initially hypothesised that not all interaction types would contribute equally to the maintenance of species coexistence. Accordingly, we observed that adding the effect of different types of interactions to the baseline of plantplant interactions (scenario 1; Fig. 3.1) did not further modify the structural niche and fitness differences (SND and SFD) that determine the triplets’ feasibility (Fig. 3.3). Analyses for modules of four species (i.e., quadruplets) rendered similar results (Fig. S3.3). For instance, incorporating the full complexity of interactions, which include the direct effects of grasshoppers on plant species plus the plantplant and plant-herbivore HOIs (scenario five), produced a slight and positive 79 3. Self-limitation and facilitation in a complex herbivore-plant network increase (0.5—1%) of SND and SFD. These results identify the structure of plantplant interactions ( α matrix) as the key element that modulates the structural mechanisms of multispecies plant coexistence and diminishes the importance of multitrophic interactions. Given the previous result, we selected the matrix of plant interactions to analyse what particular types of network configurations explain the variation observed in SND and SFD. Two particular metrics, diagonal dominance and the positive-negative ratio (PNR) usefully described SND variation (Fig. 3.4, Table S3.2). In particular, an increase in SND was observed in modules with higher PNR (stronger positive than negative values of interactions) and higher diagonal dominance (stronger strength of intraspecific interactions compared to interspecific ones). Additionally, we observed similar results in modules composed of four species, but for those, we also found significant interactive effects of both PNR and diagonal dominance on describing the observed variation in SND (Fig. S3.4 and Table S3.4). No other relevant relationships were found between SND and any other investigated metric in triplets and quadruplets. On the other hand, several network metrics and their interactions explained the observed variation of SFD across triplets and quadruplets. Regarding the PNR, we observed that stronger positive than negative interactions notably decrease SFD. Also, higher kurtosis is associated with lower SFD, which indicates that weaker interactions decreased SFD. Nevertheless, the dominance of the diagonal (ratio of intraversus interspecific effects) was not observed to have an influence on SFD (Fig. 3.5). For the case of interactions between network metrics, we observed that all two-metric combinations (except the interaction between diagonal dominance and skewness) were relevant to describe SFD variation. Finally, triple interactions (diagonal dominance—skewness—PNR and diagonal dominance—kurtosis—PNR) slightly decreased and increased SFD, respectively (Table S3.3). Very similar results were observed for modules of four species, except for the case that diagonal dominance did positively contribute to explaining SFD variation in quadruplets (Fig. S3.5, Table S3.5). Overall, these results suggest that the SFD that drive competitive dominance can be promoted by multiple network configurations of species interactions whereas fewer particular network characteristics can explain the SND that drive the opportunities for species to coexist. 80 3.4. Discussion 3.4 Discussion Ecologists have long sought to understand the role of species interaction in promoting the maintenance of biodiversity. Traditionally, this aim has been focused on systems of low diversity where estimating the strength and sign of biotic interaction that describe population dynamics was challenging but not impossible. However, the goal of recent work is to extend these efforts to highly diverse and multitrophic systems in which multiple interactions are simultaneously occurring (Levine et al., 2017; Godoy et al., 2018). It is clear that these new opportunities bring new challenges and it is meaningless to consider relevant all biotic interactions occurring between species of the same and different trophic levels, and between pairs or multiple species. Consequently, we first explored which subset of species interactions describes the spatiotemporal distribution of plant cover in our study by means of a regularisation approach. Similar methodologies have been followed in other studies working with annual plant species (Weiss-Lehman et al., 2022) and with bacteria (Coyte et al., 2021) to address this challenge. We found statistical support for all types of interactions including pairwise plant-plant and plant-grasshoppers interactions as well as HOIs, but it is also important to highlight that we found a strong pruning in the proportion of interactions retained. For instance, more than half of all potential interactions within plant species were considered to be negligible, and the proportion went up to 2.00 ± 0.89% for complex HOIs describing how the interaction between a pair of plant species is modified by the presence of a third (plant or grasshopper) species. With such information at hand, we could in principle conclude that all interaction types are relevant and what matters is maintaining a combination of simpler and more complex interactions among a subset of species. Yet, the effort we did in our study to connect this subset of interactions with the determinants of multispecies coexistence gives further insights. Specifically, we found that not all interaction types selected by the regularisation approach contributed similarly to the mechanisms determining the feasibility of the community. In particular, the plant-plant network of interactions had a predominant role in promoting the structural niche and fitness differences that drive multispecies coexistence compared to the plant-herbivore interactions and to the HOIs. The underlying explanation for this finding is that the variation in the strength and sign of interaction coefficients among plant species was larger compared to the rest of the interaction types. For instance, while for plant-plant interactions the strength varied from -0.40 to 0.45, the rest of the interactions were weaker and peaked more towards zero. Our results, therefore, agree with previous findings (Rossberg et al., 2011) and discussions (Barab ´ as, 2021) that the variation in interaction strength plays a key role in driving the feasibility of ecological systems and suggest that the wider the variation in interaction strength in a particular trophic level, the greater its importance for modulating the diversity of the system. Together with this main result, we add the novelty that different structures of interactions also promote different coexistence mechanisms. On the one hand, we observed for modules of three and four species that the opportunities for species to coexist (i.e., structural niche differences) increased when intraspecific exceeds interspecific plant competition, and when positive interactions were more 81 3. Self-limitation and facilitation in a complex herbivore-plant network prevalent than negative ones. That self-limitation processes coming from using particular soil resources or specific plant enemies are an important driver of diversity had been well acknowledged before (Chesson, 2000; Chesson and Kuang, 2008; Barab ´ as et al., 2016), but less clear is the fact that facilitation can be as important as strong intraspecific competition for multispecies coexistence (see Losapio et al. (2021) for a suggestion without including interaction strength). On the other hand, our results show that the variation in interaction strength that determines plant competition dominance (i.e., structural fitness differences) can be structured in many different ways. For instance, modules with stronger plant interactions, represented by a small kurtosis, are likely to find a stronger competitor that reduces the ability of the system to coexist. In addition, our results indicate that this small kurtosis interacts with the positive-negative ratio (PNR) suggesting that higher dominance occurs when interactions are biased towards competition but not facilitation. Taking this into an applied context, our results highlight that we can be certain in the specific ways by which we can promote or maintain diversity, that is strong self-limitation and equally important facilitation, yet competitive dominance is harder to control or predict because there are many different ways by which some plant species can be more competitive than others. Regardless of the interaction types (pairwise versus HOIs) and trophic position (plant versus grasshoppers), our regularisation approximation indicates that species are simultaneously involved in negative and positive interactions. This is an important finding for discussing how we should approximate the study of the complexity of interactions in multitrophic systems. As we introduced, the traditional approach has been to assign a priori the sign of interactions (e.g., competition and herbivory with negative effects while pollination with positive effects) (Chesson and Kuang, 2008; Bastolla et al., 2009; Rohr et al., 2014). However, these decisions have been made due to theoretical reasons and knowledge limitations, but have no empirical justifications. If we want to accurately describe a network of multispecies interactions in multitrophic systems that connect with the spatiotemporal dynamics, we advocate that we should expect positive and negative interactions across ecological communities. In other words, there are no specific types of interactions exclusively negative or positive. The most evident case in our study system comes from the layer of interactions describing the grasshopper effect of plant cover. Herbivory is generally considered to have a negative effect on plant growth, yet we found that while some grasshoppers strongly reduced plant cover in some plant species (e.g., Chorthippus dorsatus had a strong negative effect on Bromus erectus) others did the contrary (e.g., Pezotettix giornae had a strong positive effect on Galium verum). We could not test the ultimate mechanisms by which these opposing effects are occurring but according to previous studies in the same system, the sign of these effects could be related to the feeding traits of grasshoppers. It has been shown that our six grasshopper species differ in their feeding niche as a function of their incisive strength and two simple plant traits (C/N ratio and leaf dry matter content) (Deraison et al., 2015a). Another possibility is that the negative or positive herbivory effects could be related to an interactive effect between grasshopper traits and plant abundance. Bromus erectus is a dominant plant species in the system and therefore it is an abundant resource for many grasshopper species that have enough incisive strength to eat 82 3.4. Discussion these leaves with a high C/N ratio. These linkages suggest that this is a suitable system to study the role of plant and grasshopper traits in mediating the role of species interactions in the maintenance of plant diversity. The most important findings of our study rely on the regularisation approach we used to infer the structure of interactions. This approach was fitted with information of three botanical surveys and with grasshopper abundance that was held constant across the experiment. Despite the fact that this experimental effort represents a unique example in the literature, we do not certainly know whether the vast amount of data collected is enough to represent the dynamics of the system. For instance, the estimation of interaction strength and sign obtained could be different if the experiment was conducted during a longer period with more frequency of drier and wetter events. That is, we did not consider the role of environmental variation but this is something needed for future work as many studies have highlighted, for example, the importance of precipitation variability in maintaining plant diversity (Hallett et al., 2019). Another limitation of our approach is that we could study the feasibility of the plant species but not the combination of plant and grasshopper species altogether. This was due to the fact that the number of grasshoppers was held constant across the entire duration of the experiment, and therefore we could not know how plant species affect changes in grasshopper abundance across time. This information is harder to obtain from experiments conducted under natural conditions but it could be more easily done in experimental approximations like ours in which insects are kept in large cages (see Bartomeus et al. (2021) for an example with pollinators). This discrepancy between what was done and what could have been done to obtain the maximum profit from the extensive experimental effort also calls for a better dialogue between empiricists and theoreticians. Although our study presents a complex approximation to study the mechanisms based on species interactions that determine the dynamics of ecological communities, we reach a straightforward conclusion. Pairwise plant-plant interactions and the prevalence of strong plant self-limitation and facilitation are the key components of the feasibility of our system. Despite studying a plant-herbivore community, herbivores play a minor role in plant dynamics, and the plants can be seen as the substrate over which they can display their niches according to their feeding strategies. We do not know numerically whether these asymmetries in the importance of interactions across trophic levels are common in nature, as this topic has been rarely explored in the literature. Therefore, that is a clear knowledge gap that needs further exploration. Another important conclusion is that, in order to realistically address the challenge of understanding how diversity is maintained, we need to explore multispecies, multitrophic systems, and for this exploration is fundamental a combination of empirical data and statistical approaches that allow distinguishing important interactions to reproduce the spatiotemporal patterns observed. For this study, we have used state-of-the-art statistical techniques via hierarchical group lasso regularisation (Lim and Hastie, 2015) but this is an area of rapid development, which will likely benefit studies addressing the complexity of ecological communities. Taken together, our results show that simple processes of ecological interactions (self-limitation and facilitation) structure the dynamics of complex ecological communities, in our case, a 83 4. Precipitation, traits, and niche variability 90 4.1. Introduction 4.1 Introduction The concept of the niche is a central pillar in ecology. Although the definition of the term has been strongly debated over the last century, ecologists agree on the idea it is a species property that defines the range of environmental conditions a species could endure (Grinnellian niche) (Grinnell, 1917), their role within food webs (Eltonian niche) (Elton, 1927), and more generally, the multiple abiotic and biotic dimensions that allow species to develop, reproduce and persist (Hutchinsonian niche) (Hutchinson, 1957). During most of the 20th century, the definition of the niche was used to explain how biodiversity can be maintained within ecological communities. It was formulated that the biodiversity we observe is the result of evolutionary processes that have promoted enough niche differences between species to avoid competitive exclusion (Gause, 1934; MacArthur and Levins, 1967). The main prediction was that a community can accommodate as many species as limiting resources unique to each species exist. That is, if two species compete for similar resources but one is more limited by carbon (C) and the other species more limited by nitrogen (N), both could coexist without excluding each other. However, the first attempts to test these predictions showed that niche differentiation for limiting resources was not able to accurately explain the co-occurrence of multiple competing species. Some of these contradictory results between theoretical expectation and empirical findings, formulated as the paradox of the plankton, demonstrated that it is possible to observe high diversity of planktonic organisms despite being assembled in environments with few limiting nutrients, and indeed, fewer resources than the number of competing species (Hutchinson, 1961). These discrepancies stimulated a plethora of further studies highlighting that, as important as the number of limiting resources is, there are also other critical aspects of niche differentiation between species that promote coexistence such as natural enemies, rates of resource uptake and their spatial and temporal variation (Holt, 1977; Chesson and Warner, 1981; Chesson, 1985; Tilman, 1990, 1993). The advances in mechanistically understanding the multiple abiotic and biotic sources of niche differences reinforced the idea that the niche has a multidimensional nature (Chase and Leibold, 2003). Yet, these advancements were still not useful to predict if species could coexist, or even if they coexist, whether species are rare or abundant or whether their abundances could change with time. The link between the species’ demography with their niche differences was done by P. Chesson during the decades of the 80s, and 90s. His contributions, summarised in Chesson (2000), posit that because knowing all sources of niche differences between species is almost impossible, we can instead focus on the demographic signature of niche differences. This demographic signature allows for maintaining the diversity of species by limiting their abundance when they become common but also by buffering them from extinction when they become rare. Nevertheless, this framework, currently named modern coexistence theory (MCT), brought another lesson: niche differences alone do not determine whether species can coexist, they only stabilise the dynamics of interacting species. There is another source of species differences. These other differences, named fitness differences, arise from the combination of the intrinsic ability of species to produce 91 4. Precipitation, traits, and niche variability offspring and their tolerance to competition, promote competitive dominance, and in the absence of niche differences indicate which species is the superior competitor (Godoy and Levine, 2014). With both species differences acting simultaneously, MCT predicts that species can coexist when niche differences “overcome” fitness differences (Chesson, 2000; Barab´ as et al., 2018). Applications of MCT during the last two decades have triggered an unprecedented advancement of our mechanistic understanding of the role of competition (Allesina and Levine, 2011; Hart et al., 2016, 2018; Adler et al., 2018), shared natural enemies (Chesson and Kuang, 2008; Mora et al., 2011), mutualisms (Lanuza et al., 2018; Johnson et al., 2022), and environmental variability (Mat ´ ıas et al., 2018; Hallett et al., 2019) on maintaining biodiversity through promoting niche differences. Yet, mathematical and methodological tools associated with MCT cannot measure the niche of a species but only differences between pairs of species. Fortunately, recent progress by Spaak and De Laender (2020) has provided an intuitive and broadly applicable definition of the niche and the fitness that can be evaluated at the species level. These new advances allow us now to answer questions such as which is the position of one species with respect to the rest of the community (Spaak and De Laender, 2020), but rather more importantly, it provides the tools for asking, for the first time, whether the niche is an intrinsic property of the species. If that is the case, then niches could be differentiated by analysing differences among species and differences also in their functional characteristics that determine the constraints of the niche. Alternatively, it could be also possible that the niche is an ecological property more variable than previously thought, and the variation in niches we find within ecological communities can be better explained by environmental variation or by species abundances rather than by their identity. This is for instance, the case of the species’ fitness, that has been shown to operate as a mutable property determined by environmental drivers (Mat ´ ıas et al., 2018; Wainwright et al., 2019; Van Dyke et al., 2022) and species abundances (Hart et al., 2016), which result in competitive rankings among species to be susceptible to change across space and time (Stouffer et al., 2021). Here, we challenge the central idea in ecology that the niche is a species’ property. To do so, we applied the most recent, species-level definitions of Spaak and De Laender (2020) to quantify the niche and the fitness of 18 annual plant species in Mediterranean grasslands for five consecutive years. This temporal variation allows us to test whether the niche position, and therefore niche differences between species were consistently maintained among species across time, or conversely, the niche position of species we observe does not carry taxonomic signals. Following current theory, we would expect greater variability in species fitness than species niche, with the fitness varying greatly between years depending on environmental conditions and the niche behaving in a more conservative way, showing less variation and being mainly explained by the traits of species that allow them to thrive under a limited range of environmental conditions. 92 4.2. Materials and methods 4.2 Materials and methods 4.2.1 Study system and fieldwork Our experimental site (Caracoles Ranch) is located within Do ˜ nana National Park in southwest Spain (37° 04’ N, 6° 18’ W). The area has a Mediterranean climate with mild winters and an average annual precipitation of 550-570 mm over the past 50 years. The physiognomy of the vegetation is primarily grasslands of annual grassland species, with no perennial species present. Across the experimental site, there is a slight slope of 0.16%, which causes vernal pools to form at the lower border of the ranch during the winter and spring months (November-January and March-May) while the upper parts of the ranch do not flood, except during exceptionally wet seasons. During the course of our study, an extreme flooding event occurred during the growing season of 2018. There is also a strong soil salinity-humidity gradient that is structured along this topographic gradient. In September 2014, we established nine plots of 8.5m x 8.5m along an area measuring 1km × 200m. The plots were divided into three groups, with three plots located at the upper part of the topographic gradient, three in the middle, and the remaining three at the lower part. The average distance between these three groups of plots was 300m and the average distance between plots within each group was 30m (with a minimum distance of 20m). Each plot was further divided into 36 smaller plots of 1m × 1m, with 0.5m aisles in between for access during measurements (for a total of 324 subplots). This layout was implemented to facilitate the parameterisation of population models with which we could estimate the species niche and fitness following Spaak and De Laender (2020) definitions. These models include two types of parameters. On one hand, the intrinsic fecundity component that describes the average intrinsic ability of species to produce seeds, which is in turn discounted by the seed amount that is lost due to germination and due to seed mortality in the soil (Godoy and Levine, 2014). On the other hand, the effect of intra- and interspecific pairwise interactions that describe how this intrinsic ability of seed production is reduced by an increase of conspecific and heterospecific densities. In order to parameterise these models, the main focus of the field observations involved measuring, for each focal individual, the production of per germinant viable seeds as a function of the number and identity of neighbours within a radius of 7.5 cm, including individuals of the same species. This radius is a standard distance used in previous studies to measure competitive interactions among annual plant species (Levine and HilleRisLambers, 2009; Mayfield and Stouffer, 2017) and has been validated to capture the outcome of competition interactions at larger scales (1m 2 ) under locally homogeneous environmental conditions (Godoy and Levine, 2014). From November 2014 to September 2019, we sampled 18 species present in the study area each year (Table S4.1). There were two types of species, those considered core to be a species found every year (8 species) and non-core to be a species that was missing at our plots but not at our experimental site for at least one year (10 species) (Table S4.1). We sampled one individual per subplot for widespread species and several individuals per subplot when species were rare (max. 324 individuals/species). This sampling design ensured that all species were balanced in terms of the 93 4. Precipitation, traits, and niche variability number of observations and that we captured the full range of observed spatial interactions among species across the study area. Furthermore, we obtained independent estimates of seed survival in the soil and seed germination rates in 2016 (see Lanuza et al. (2018) for details on obtaining these rates). 4.2.2 Estimating species interactions We obtained pairwise per-capita interaction coefficients, both intra- and interspecific, from the empirical data in each year of sampling. In particular, we modelled the observed viable seed production per individual (yield of species i or Yi ) as a function of the identity and abundance of neighbouring individuals, following a previous study in the same study area (Garc ´ ıa-Callejas et al., 2021). The per-capita effect between each pair of species was estimated by assuming an underlying Ricker model of population dynamics (Ricker, 1954; Mayfield and Stouffer, 2017) of the form: Yi=λie−Pn j=1(αij Nj)(4.1) where λi is the number of seeds produced by species i in the absence of interactions, αij is the per-capita effect of species j over species i (which can be positive or negative, thus allowing both competitive and facilitative effects), and Nj is the number of individuals of species j within 7.5 cm of the focal individual. We implemented this model in a Bayesian framework, via a multilevel model with a negative binomial distribution and allowing the intercept and slope of the relationships to vary across years but not across space by including the year as a random effect. Thus, the αij values in eqn. 4.1 vary across years but are homogeneous across the whole spatial scale for the whole study area. We used the means from the obtained posterior distributions of αij values as the coefficients of our interaction matrices. For fitting the model, we used non-informative priors with MCMC settings of 5,000 iterations (of which 2,500 were warm-up) and 6 chains. We implemented the model using the “brms” R package (B¨ urkner, 2017). 4.2.3 Calculation of niche and fitness With the empirical estimations of species vital rates and interaction coefficients previously described, we followed recent definitions by Spaak and De Laender (2020) to quantify the niche and fitness at the species level following previous intuitions (Adler et al., 2007). These estimations are defined within the framework of modern coexistence theory (MCT) and are also framed within invasion analyses (Grainger et al., 2019), but they differ from the traditional approach, among other things, because they are no longer calculated as a difference between pairs of species (Chesson, 2000). The advantage here is that by defining the niche at the species level, we can explore whether such estimates vary more between species than between years and whether they are predicted by species traits (indicating that it is a species property) or vary more across years than across species (indicating that it is not a species property). 94 4.2. Materials and methods The definitions of Spaak and De Laender (2020) assume that each species has a stable monoculture equilibrium and that the invasion growth rate correctly predicts coexistence. That is, for coexistence to occur, all species need to have positive invasion growth rates, defined as the growth rate of a species when it is reduced to density near to zero and the other species are at their monoculture equilibrium. The calculation of niche differences for a multispecies community is made as follows: Nichei= Invasion growth rate z }| { fi(0,N−i,∗)− No-niche growth rate z }| { fiX j,i cijN−i,∗ j,0 fi(0,0) | {z } Intrinsic growth rate −fiX j,i cijN−i,∗ j,0 | {z } No-niche growth rate ,(4.2) where fi is any function that describes the per capita growth rate of species i ; the term N −i,∗ is the vector of equilibrium densities in the community when species i is absent; cij denotes a conversion factor that can be calculated empirically to account for the difference in resource consumption of individuals from the different species (e.g., the difference in resource use between a tree and a herb), translating the density of a species j into the density of species i ; and 0denotes the absence of all species other than i. The framework introduces a term called no-niche growth rate of species i , which corresponds to the growth rate of species i if there was no niche differentiation, i.e., the growth rate of i at the converted monoculture density of its competitor species j (Spaak and De Laender, 2020). The niche is calculated as a comparison between the growth rate of species i when only interspecific interactions are present (invasion growth rate) and the growth rate when only intraspecific interactions are present (no-niche growth rate), divided by the comparison of the growth rate when the density of species iis close to zero (intrinsic growth rate) and the density of i is at the converted equilibrium density of species j . In other words, the niche measures the strength of the frequency dependence divided by the strength of the density dependence of species i (Adler et al., 2007). With this in mind, the fitness of species i can be calculated for the scenario in which all species occupy the same niche (Adler et al., 2007; Barab ´ as et al., 2018), thus measuring how well species i grows in the absence of frequency dependence, i.e., when there are only intraspecific interactions (no-niche) against its intrinsic growth rate: Fitnessi=fiPj,icijN−i,∗ j,0 fi(0,0).(4.3) 95 4. Precipitation, traits, and niche variability Both (eqn. 4.2 and 4.3) show that species abundances at equilibrium (N∗) are necessary to define the species niche and fitness following Spaak and De Laender (2020) definitions. Therefore we computed these equilibrium abundances by parameterising a population model with estimates of λi and αij obtained from eqn. 4.1. The population model, which was a discrete-time logistic model, followed the same functional form (Ricker) as eqn. 4.1. 4.2.4 Functional traits and climatic variables For each species, we measured at the beginning of the experiment (2015) during the peak of biomass of each individual species (seasonal development ranged from early February to late September), a series of plant functional traits related to particular organs such as leaves (leaf area, and specific leaf area) and roots (root diameter, specific root length, root tissue mass density, and specific root area), as well as whole-plant traits (i.e., height, canopy shape, leaf area index, C:N ratio, and C13 and N15 isotopes) (Table 4.1). These traits were selected for their recognised or assumed utility as response traits to abiotic conditions and their potential implication in plant competitive dynamics (P ´ erez-Ramos et al., 2019). Plant sampling and trait measurements followed standard protocols recommended by P ´ erez-Harguindeguy et al. (2016). Briefly, 30-50 individuals per species across their range of distribution were randomly selected to measure height and canopy shape, and leaf-level traits. Leaf size was quantified using an image analysis program (Image Pro-plus 4.5; Media Cybernetic Inc., Rockville, MD, USA). The four morphological root traits were measured in 3 to 7 individuals per species and treatment. In the field, we separated root systems from individuals of different species carefully. In the laboratory, we washed root samples to remove soil and then a representative subsample of fresh fine roots (less than 2mm in diameter) was scanned at 1200 dpi. The digital images were used to determine the length, area, volume, and mean diameter of roots using specific image-analysis software (WinRHIZO version 2003b, Regent Instruments Inc., Quebec, Canada). The root material harvested was immediately weighed, oven-dried at 60°C for 48 hours and then re-weighed. Four chemical traits (leaf N and C content, as well as the bulk isotopic composition of these elements) were determined on three composed samples per species and treatment obtained by pooling the leaves previously used for above-ground traits, using a CHN elemental analyser coupled to an isotope mass spectrometer at the laboratory facilities of Do ˜ nana Biological Station (EBD-CSIC). We also obtained precipitation information, both annual and seasonal (winter, spring, summer, and fall) in mm (l/m 2 ) for each year from the nearest weather station (Aznalcazar, 10km far apart) that belongs to the Autonomous Government of Andalusia. 4.2.5 Statistical analyses To test whether observed temporal variation in species’ niche and fitness was greater between years or between species, we performed pairwise t-tests (a t-test between each year) using BH as the adjustment method to minimise true positives 96 4.2. Materials and methods (Benjamini and Hochberg, 1995). Once we had understood the source of niche and fitness variation, we tried to link this variation to properties related to species identities (functional traits) and to environmental heterogeneity (annual and seasonal precipitation). For linking niche and fitness with species functional traits, we fitted a permutational multivariate analysis of variance (PERMANOVA) with 999 permutations and euclidean method (Anderson, 2014); the results are shown in Table 4.1. Finally, to link niche and fitness with the precipitation, we applied linear (e.g., niche ∼ rain) and quadratic (e.g., niche rain 2 + rain) generalised linear models (GLM) for each variable (total annual rain as well as total rain during winter, spring, summer, and fall) and then selected the best model based on Akaike information criterion (AIC) (Bozdogan, 1987). All computational analyses were performed in R version 4.2.1 (R Core Team, 2022). 97 4. Precipitation, traits, and niche variability 4.3 Results Figure 4.1: Relationship between species’ niche and fitness for the 18 studied plant species across five years of field observations (2015-2019). Ellipses drawn assuming a multivariate t-distribution with 95% confidence. A species is considered core when it is present in all sampled years (see Table S4.1). Differences between years following a pairwise t-test can be found in Table S4.2 for fitness differences and in Table S4.3 for niche differences. We estimated the niche and fitness for our 18 plant species (6 of which are core to the community, i.e., are found in all five sampled years; see Table S1) thanks to coupling a strong empirical effort with recent theoretical advances that define these characteristics at the species level. Contrary to the definition of the niche in any of its potential meanings, we observed that the niche is not a species property. That is, the position of the niche of a given species greatly varies across years (Figs. 4.1 and 4.2), and the relative differences between pairs of species were not temporally maintained (Fig. 4.1). Even more interesting is that we observed a directional change in the species’ niche as well as their fitness from one year to another. This directional change indicates the critical role of the abiotic environment in defining the species niche, and in allowing species to explore different combinations of niche and fitness across years due to associated temporal variation in precipitation. We also observed a great variation in species niche and fitness with almost all years rendering different results compared to the other four. For instance, only 98 4.3. Results Figure 4.2: Annual changes in species’ fitness (a) and niche (c), and (b) relationship between them. A species is considered core when it is present in all sampled years (see Table S4.1). Lines across years show changes in niche and fitness ranking for core species. Differences between years following a pairwise t-test can be found in Table S4.2 for fitness differences and in Table S4.3 for niche differences. in 2018 did species show similar fitness values to 2015 and 2016 (Fig. 4.2b and Table S4.2). Moreover, we did not find two consecutive years in which the average fitness differentiation remained unchanged. In that sense, it is worth highlighting the notable differences between 2017 and 2019, with average fitness differences of -0.08 ± 0.07 and 0.13 ± 0.05, respectively. For the case of the species’ niche, we found very similar patterns of variation (Fig. 4.2c and Table S4.3). Except for 2017, which showed similar niche values as 2016, we observed strong significant differences between the niches over the years. Particularly in 2018, the year in which our field site went through a big flood, the variation in species niche was the largest (2.60 ±0.68) compared to, for instance, 2015 (1.49 ±0.12). As a consequence of the high variability observed for our system, we could neither assign with confidence a hierarchy of competitive dominance (Fig. 4.2b) nor 99 4. Precipitation, traits, and niche variability Associated publication • Granjel RR, Spaak JW, Buche L, Garc ´ ıa-Callejas D, Arroyo J, Godoy O. Precipitation, not species functional traits, explains the interannual variability in plant niche. In preparation. Open science Data and code available through www.github.com/Granjel/Caracoles. Author contributions O.G. conceptualised the research. O.G., R.R.G., and L.B. performed fieldwork. D.G.-C. parameterised the population model. J.W.S. quantified niche and fitness. R.R.G. analysed and designed the visualisation of results. R.R.G. and O.G. wrote the original manuscript with input from all authors. Acknowledgements We thank the Radical Community Ecology group for discussions and for contributing to data acquisition at Caracoles Ranch. Additional funding O.G. acknowledges the postdoctoral financial support provided by the Spanish Ministry of Economy and Competitiveness (MINECO) and by the European Social Fund through the Ram´ on y Cajal Program (RYC-2017-23666). Supplementary material Supplementary figures and tables related to this chapter are available in Appendix 4. 106 Discussion 107 Discussion 108 The rapid pace of anthropogenic global change is dramatically altering the natural world. Shifts in environmental conditions are leading to modifications of biotic interactions such as competition, mutualism, or predation (Tylianakis et al., 2008). In turn, the balance of such interactions is responsible for determining whether species are able to persist in a community or, on the contrary, they go extinct. Therefore, understanding how environmental changes and ecological interactions are shaping the mechanisms that allow species to coexist is nowadays more pressing than ever. Throughout the chapters of this thesis, I have raised fundamental questions in community ecology and approached them from an integrative perspective that combines the most recent theoretical advances in the field of ecology with extensive empirical data. The approach used here is unmatched, as this thesis comprises the first effort towards understanding how different sources of abiotic and biotic heterogeneity affect the mechanisms of coexistence in multispecies, rather than pairwise, ecological communities. In the following lines, I discuss the main findings and developments of this thesis along with its main limitations, which open paths for future work. Finally, I present a list of the main conclusions reached as a result of this work. In this thesis, we explore the effects for species coexistence of two primary drivers of global change: nitrogen enrichment and changes in land use intensity (Tylianakis et al., 2008). On the one hand, nitrogen (N) enrichment is currently challenging our understanding of a thoroughly studied topic in ecology such as resource competition. Much of the evidence regarding N enrichment in natural ecosystems indicates that predicted increases in the availability of this nutrient are related to shifts in species composition and biodiversity losses (Stevens et al., 2018; Midolo et al., 2019). This is suggested to happen through multiple mechanisms, such as an increase in fast-growing species leading to asymmetries in the competition for light (DeMalach et al., 2017) or a reduction in niche dimensionality that increases interspecific interactions (Harpole et al., 2016). However, our results indicate that an increment of N may not always lead to biodiversity loss (Chapter 1). We showed that the experimental addition of N to perennial plant communities did not predict losses of biodiversity. Instead, it promoted the niche differences that increase the chance of persistence and decreased the fitness differences that arise from competitive abilities and potentially lead to competitive exclusion. In light of all previous research, we claim that, although this does not indicate the existence of beneficial effects of N enrichment for natural ecosystems, higher N levels have the potential to increase niche differentiation and decrease fitness differentiation in certain years, which in turn may buffer its negative impacts on biodiversity. On the other hand, land use is a complex concept comprehending different variables such as fertilisation, mowing, or grazing. Increasing land use intensity (LUI) has apparent negative consequences for biodiversity, including declines in plant diversity and compositional changes (Rockstr ¨ om et al., 2009; Newbold et al., 2015; Allan et al., 2014; Crawley et al., 2005; Pakeman, 2011; Dormann et al., 2007); however, we have little knowledge of the mechanisms behind this pattern. Based on the existing literature on the topic (e.g. Flynn et al., 2009; Lalibert ´ e and Tylianakis, 2012; Mayfield et al., 2010), we have identified three possible 109 Discussion ways by which increases in land use intensity could hinder species persistence: by increasing competitive inequalities and thus fitness differences, by reducing the axes of niche differentiation, or by limiting the capacity of indirect interactions (discussed below) for promoting biodiversity. What we found when we analysed 150 grasslands under a gradient of LUI (Bl ¨ uthgen et al., 2012) across eight years (chapter 2) is that, although fitness differences were high and restricted coexistence, they did not change in response to the gradient of LUI. Niche differences, however, show an unprecedented nonlinear relationship with LUI that implies higher chances of coexistence in lower and intermediate values of LUI than in higher values. Overall, we found relatively high niche differences across the LUI gradient. We are not able to elucidate the specific components of our LUI index that triggered such high niche differentiation, but associations with organisms that do not suffer from land-use intensification such as mycorrhizae (Hart et al., 2003), the effects of pathogens, or trade-offs with different nutrients (Harpole et al., 2016) could be plausible explanations worth exploring. Moreover, we observed high and widespread interspecific facilitation that, as pointed out in recent studies by Adler et al. (2018) and Losapio et al. (2021), could be one of the main factors behind biodiversity maintenance across environments, even in anthropised ecological systems. Within this thesis, we also explored the role of natural enemies such as foliar fungal pathogens and herbivores for species coexistence. Fungal pathogens were studied within the same experiment as N enrichment (chapter 1) because fungal pathogens can benefit from either higher N concentration in leaves (Dordas, 2008) or a shift towards fast-growing plant species (Cappelli et al., 2020) caused by N enrichment. Pathogens, as natural enemies, are thought to promote biodiversity via stabilising and equalising mechanisms, i.e., both by suppressing increasingly common species and reducing their competitive ability (Petermann et al., 2008; Mordecai, 2011; Bagchi et al., 2014; Bever et al., 2015). However, the effects of pathogens may be complex and hard to predict because changes in pathogen abundance and composition could benefit a particular combination of specialist and generalist species, potentially resulting in neutral (Spear and Mordecai, 2018; Uricchio et al., 2019) or even detrimental effects for biodiversity maintenance. Our results indicated that the application of a fungicide treatment to suppress foliar fungal pathogens caused similar effects to those of N enrichment. In other words, plant communities with natural, unmodified pathogens showed lower niche and higher fitness differences than those plant communities with suppressed pathogens. Although we lack a detailed understanding of the pathogen community to be able to decompose the specifics by which the fungicide treatment promoted the determinants of plant coexistence, one possible explanation that we mention is that the natural pathogen community may have been composed of fungal species with overall neutral effects for coexistence (e.g., not engaging in density dependence mechanisms). Thus, applying the fungicide may have selected generalist species prone to attack the most common species in the community (see Cappelli et al. (2020)). As well as pathogens, insect herbivores are expected to promote coexistence through negative density dependence (Ishii and Crawley, 2011; Heard and Sax, 2013; Forrister et al., 2019; Descombes et al., 2020). However, their effects are also 110 susceptible to change depending on the composition of the herbivore community and their feeding preferences (Kempel et al., 2015). To understand this, we analysed an empirical network of perennial plants and herbivore grasshoppers (Chapter 3). For this network, we estimated the direct interactions between plants and the effects of grasshoppers on plants, but also we estimated the higher-order interactions (HOIs; modifications of a direct interaction by a third species) caused by plants and by grasshoppers; such HOIs are thought to stabilise competitive dynamics and explain unexplained variability in ecological communities (Grilli et al., 2017). What we found is that most of the variation in the plant’s niche and fitness differences was explained solely by the direct interactions between plants and that neither the direct effect of herbivores nor the HOIs played a major role in species coexistence in our system. This result may be influenced by unknown determinants intrinsic to our system, but there is little doubt that this result contrasts with most previous theoretical and empirical findings and deserves further exploration. A key decision for this study was to apply state-of-the-art statistical techniques via hierarchical group lasso regularisation (Lim and Hastie, 2015), for two reasons. One, because it helped us select the important interactions explaining the dynamics of plant species within our complex plant-herbivores network, instead of trying to quantify all interactions independently of their importance. Second, this approach also allowed us to relax our assumptions regarding the sing of effects of herbivores over plant species. As we argued, the sign of interactions had traditionally been assigned aprioristically, with herbivores being typically negative and pollinators being typically positive for plant growth (Chesson and Kuang, 2008; Bastolla et al., 2009; Rohr et al., 2014). However, we observed both positive and negative effects across ecological interactions in our network, including the effects of herbivores on plants. Thus, for interactions such as those between insect herbivores and plants, where there may be hidden mechanisms that we cannot decipher, we advocate for letting ecological models determine the sign of the effect of interactions on species dynamics. Within the same system, we explored a key question in community ecology: how does the configuration of species interactions affect the determinants of species coexistence? To do so, we used different metrics related to the distribution of species interactions within modules, to the dominance of intraspecific over interspecific interactions, and to the prevalence of positive over negative interactions. The stabilising effect of dominant intraspecific interactions within modules was already well known within modern coexistence theory (Chesson, 2000; Barab ´ as et al., 2016), and we corroborated it in our system. However, we found that the prevalence of stronger positive than negative interactions also promoted the axis of niche differentiation between perennial plants, which had been suggested (Losapio et al., 2021) but not explicitly shown in empirical weighted networks. As opposed to the niche differences, the axis of fitness differentiation could be explained by many more metrics, alone or interactively. For example, higher dominance and high kurtosis occurred mainly when competitive effects were prevalent (as determined by a low ratio between positive and negative interactions), representing increased fitness differences and higher chances of encountering a competitive dominant species destabilising the system. All of 111 Discussion this suggests that niche differentiation in ecological systems may be easier to control than fitness differentiation because fewer predictors are able to capture its variability. Applying a structural stability approach to species coexistence in Chapters 1, 2, and 3 has allowed us to explore the effects of all the topics discussed above (N enrichment, land use intensity, pathogens, herbivores, and interaction configurations) on the determinants of multispecies coexistence. One principal result of this thesis is that the number of interacting species determines the effects of different sources of abiotic and biotic heterogeneity over the niche and fitness differences that determine species coexistence. Specifically, N enrichment and pathogen suppression effects on plant coexistence (Chapter 1) differed between pairs and multispecies assemblages (three to six); changes in both niche and fitness differences were more pronounced with the gradient of LUI (Chapter 2) in three-species assemblages compared to pairs of species; herbivory and HOIs (Chapter 3) had even less relevance explaining the variation of niche and fitness differences in four-species combinations compared to three-species ones; and we found interactive effects between positive interactions and intraspecific interactions explaining niche differences in four-species combinations as opposed to three-species ones (Chapter 3). The exploration of multispecies mechanisms requires further efforts to decipher the specific causes of why different numbers of interacting species show different responses to abiotic and biotic heterogeneity. However, with only a few available studies accounting for multispecies coexistence mechanisms (Petry et al., 2018; Tabi et al., 2020; Bartomeus et al., 2021; Deng et al., 2021; Garc ´ ıa-Callejas et al., 2021; Johnson et al., 2022; Arroyo-Correa et al., 2023), the empirical results presented here are very novel and suggest that we need to explicitly incorporate the number of interacting species into coexistence studies if we want to accurately predict the effects of global change drivers and biotic interactions on biodiversity maintenance. Additionally, another main benefit of applying a structural approach is that we are able to quantify the contribution of indirect effects, which arise in multispecies assemblages, to the determinants of species coexistence (Saavedra et al., 2017). One of the best-studied indirect effects in ecology is intransitivity, represented by the popular game of rock-paper-scissors dynamics in which rock beats scissors, scissors beat paper, and paper beats rock, stabilising an otherwise unstable community (Gallien et al., 2017; Soliveres et al., 2018). These effects are evaluated in Chapter 1 and Chapter 2 of this thesis. We found that indirect effects created new opportunities for species persistence under N enrichment and pathogen suppression (Chapter 1), but these opportunities were relevant at simpler species combinations (three and four species) and they lost importance as more species were involved (combinations of five and six species). Moreover, we did not predict changes in the importance of indirect interactions with the gradient of LUI (Chapter 2), which aligns well with previous results indicating a minor presence and role of intransitive dynamics in annual plants (Godoy et al., 2017; Mat ´ ıas et al., 2018) but contrasts with previous results in the study system (Soliveres et al., 2018). This may be explained by a neutralising effect of alternative indirect effects such as indirect facilitation, which occurs when a species strongly suppresses a shared competitor with a third species (Levine, 1999; Flory and Bauer, 2014). Im- 112 portantly, the metrics used were useful for describing the contribution of indirect interactions to coexistence, described our empirical results effectively, and have the potential to be used in prospective studies. In the three first chapters, we embraced the latest developments in coexistence theory, structural stability, which uses structural analogues of niche and fitness differences to understand multispecies persistence (Saavedra et al., 2017; Song et al., 2018). However, we did so without questioning the nature of such important underlying concepts. In particular, the niche has been debated for decades, and today it is thought to be a multidimensional species property accounting for all abiotic and biotic heterogeneity that makes a species thrive (Chase and Leibold, 2003). But we argue that, despite being considered a species property, the frameworks that use the niche to explain crucial ecological features such as species diversity (Chesson, 2000; Adler et al., 2007; Saavedra et al., 2017) are not able to quantify it at the species level, and instead calculate differences between pairs of species or even multispecies assemblages. Thus, we wondered whether the niche actually functions as a species property or it rather behaves as a characteristic that varies with environmental conditions. In the quest to apply the most recent theoretical developments to answer central ecological questions, we applied new definitions in modern coexistence theory that allow us to quantify species’ niches and fitness at the species level (Spaak and De Laender, 2020). In our study (Chapter 4), we saw that the plant niche has a strong interannual variability that relates neither to species identity nor to species functional traits. Instead, this temporal variability can be explained by an environmental variable such as seasonal precipitation, reflecting how dry or wet each season is. Although the results of this study may be limited to short-living organisms such as annual plants, which renew their population each year and thus may be more susceptible to changes in niche and fitness, our results are robust in terms of niche variability and its lack of relation to species-level characteristics. Despite the solid combination of theory and empirical data, the results reached in this thesis come with some limitations that are worth discussing. First, in this thesis, we have not intensively analysed the effect of spatiotemporal heterogeneity. In some studies, we have multiple years of data or different locations; for example, in Chapter 1 we have two years of experimental data and in Chapter 2 we have eight years and 150 grasslands. However, such variability is not accounted for and is rather used to improve the estimation of intrinsic growth rates and interaction coefficients between species. Multiple studies show that spatiotemporal heterogeneity is key in determining species persistence in the long term through various mechanisms (e.g., Chesson, 1985; Chesson et al., 2014; Shoemaker and Melbourne, 2016; Garc ´ ıa-Callejas et al., 2021), so explicitly accounting for such heterogeneity would be key in following empirical efforts to understand the mechanisms of multispecies coexistence. Second, although this thesis contributes to understanding how multiple key ecological factors such as natural enemies or resources drive coexistence, it does not include some fine-scale details that could be determinant to acquiring a more mechanistic understanding of these processes. For example, in Chapter 1 we suppressed fungal pathogens using a fungicide but we did not further analyse the changes in the pathogenic community; the resulting composition of such community, as we have previously indicated, could 113 Discussion be critical to linking density-dependent mechanisms that stabilise coexistence to individual types of generalist or specialist pathogens and their associations with plants species (Mordecai, 2011, 2015; Cappelli et al., 2020). Another example is Chapter 2, where we analysed the effects of land use intensity as a single index that summarised the different factors involved (Bl ¨ uthgen et al., 2012); as we have indicated, our research could not investigate the individual effects of the components of the land use index (mowing, grazing, and fertilisation) because of their high correlation (DeMalach et al., 2017), but we would benefit from deciphering the joint contribution of these different components in our system, but also of their interactions with relevant biotic drivers such as mycorrhizae (Gossner et al., 2016; Hart et al., 2003) or plant pathogens Harpole et al.. The results of this thesis would also have benefited from including more information about the different processes occurring at different development stages of a plant’s life cycle. In our studies with perennial (Chapters 1, 2, and 3) and annual (Chapter 4) plants, we have focused on the interactions between adult individuals, ad we have not included competition, predation, or other relevant factors influencing the growth and survival of seeds and seedlings, although these could hold the potential to impact plant dynamics as much as phenomena occurring at later stages (e.g., Hanley and Sykes, 2009; Van Couwenberghe et al., 2013; Petry et al., 2018). Future research opportunities are varied and extensive. The methodology used in this thesis provides the tools needed to advance our knowledge of how biodiversity is maintained. An intuitive path to follow would be to investigate ecological systems including multiple types of trophic layers, such as predator, floral visitors, or mycorrhizal associations. Moreover, the key to acquiring a global knowledge of species persistence in the system would be to quantify species interactions and growth rates not only for one trophic layer but for all, as described by (Godoy et al., 2018). The experiments and observations needed to do this are often challenging or logistically unfeasible, but some approaches have already started with simplified systems kept in experimental enclosures, such as (Bartomeus et al., 2021). Another exciting future path would be not only to explore how biotic interactions stabilise ecological communities in variable environments but also to try and identify those key biotic interactions and coexistence mechanisms that optimise the function of the system (Godoy et al., 2020), which would be critical to acquiring a global understanding of the relationship between biodiversity and ecosystem functioning. The knowledge gained from and tools applied in these perspectives can be applied in interesting ways, such as helping ecologists enhance conservation and restoration efforts (Aoyama et al., 2022). The ultimate goal of this thesis is to contribute to the scientific understanding of the mechanisms maintaining biodiversity in complex and diverse communities, serving as a valuable resource for future research efforts in the field. 114 Conclusions • The use of a structural stability approach allows us to accurately predict the direct effects of global change drivers and biotic interactions on biodiversity maintenance and also to quantify the contribution of indirect effects emerging in multispecies assemblages. Future studies attempting to use empirical data to understand multispecies coexistence and ecosystem stability would benefit greatly from its application. —Chapters 1, 2, and 3. • Nitrogen enrichment is considered a global change driver related to species losses from terrestrial ecosystems worldwide, yet in our study it has the ability to increase the niche differences that promote plant coexistence and also decrease the fitness differences that lead to competitive exclusion, potentially buffering its negative impacts on biodiversity. —Chapter 1. • Suppressing the community of foliar fungal pathogens through the application of a fungicide is expected to reduce plant coexistence through the elimination of natural enemies. However, we observed similar effects as nitrogen enrichment in reducing niche and increasing fitness differences among multispecies assemblages. This unexpected result could only be obtained by the combination of theory with detailed empirical experiments, suggesting either that the natural pathogen community has neutral effects on coexistence or that the fungicide application may have selected generalist species that attack the most common plant species. —Chapter 1. • The effects of nitrogen enrichment and foliar fungal pathogens over the structural coexistence determinants depended on the number of interacting plant species, only emerging in multispecies assemblages and not in species pairs. Also, the higher the number of interacting species, the fewer species combinations are predicted to coexist. Although nitrogen enrichment and pathogen suppression do not affect coexistence predictions, they trigger changes in community composition. All of this highlights the importance of understanding coexistence beyond pairwise interactions. —Chapter 1. • Land use intensity, a complex global change driver encompassing different intensive farming practices, does not change fitness differences between plants but it does change niche differences in a nonlinear way, implying higher chances of coexistence in lower and intermediate values of land use intensity. —Chapter 2. • Interspecific facilitative interactions are strong and widespread across the gradient of land use intensity observed. This facilitation seems to play a key role in maintaining biodiversity. —Chapter 2. • Indirect effects emerging in multispecies systems promote coexistence but their effect is only noticeable when a few species interact; this indicates that, in more diverse communities, direct effects are the main force shaping 115 Appendix 1 — Chapter 1 Supplementary Figures Figure S1.1: Relationship between biomass (g) and cover (%) for all 8 species together and individually, in (a) 2017 and (b) 2018. The R2coefficient is shown for each linear model. 122 Figure S1.2: Change in biomass (g) in 2017 and 2018 for the 8 studied plant species: (To) Taraxacum officinale, (Cb) Crepis biennis, (Ra) Rumex acetosa, (Dg) Dactylis glomerata, (Ao) Anthoxanthum odoratum, (Cj) Centaurea jacea, (Sp) Salvia pratensis, and (Pm) Plantago media. The dashed horizontal line denotes the limit between positive and negative change in biomass. Note that Rumex acetosa in 2018 is the only species with a negative change in biomass. 123 Appendix 1 — Chapter 1 Figure S1.3: Structural niche differences (SND) for all combinations of two-six species, in the four experimental treatments (control, nitrogen, fungicide, and combined), for the years (a) 2017 and (b) 2018. Asterisks above boxes indicate significant differences (p < 0.05) between the corresponding treatment and the control (Table S1.9—S1.10). Note that the axis scales differ between panels because results can be compared across treatments within each species combination but not between species combinations (see main text). 124 Figure S1.4: Structural niche differences (SND) for all combinations of two to six species, in the four experimental treatments (control, nitrogen, fungicide, and combined), calculated using the coefficients estimated for the (a) lower and (b) upper limits of the 95% confidence interval. Asterisks above boxes indicate significant differences (p < 0.05) between the corresponding treatment and the control (Table S1.13 and Table S1.15). Note that the axis scales differ between panels because results can be compared across treatments within each species combination but not between species combinations (see main text). 125 Appendix 1 — Chapter 1 Figure S1.5: Structural fitness differences (SFD) for all combinations of two-six species, in the four experimental treatments (control, nitrogen, fungicide, and combined), for the years (a) 2017 and (b) 2018. Asterisks above boxes indicate significant differences (p < 0.05) between the corresponding treatment and the control (Table S1.11—S1.12). Note that the axis scales differ between panels because results can be compared across treatments within each species combination but not between species combinations (see main text). 126 Figure S1.6: Structural fitness differences (SFD) for all combinations of two to six species, in the four experimental treatments (control, nitrogen, fungicide, and combined), calculated using the coefficients estimated for the (a) lower and (b) upper limits of the 95% confidence interval. Asterisks above boxes indicate significant differences (p < 0.05) between the corresponding treatment and the control (Table S1.14 and Table S1.16). Note that the axis scales differ between panels because results can be compared across treatments within each species combination but not between species combinations (see main text). 127 Appendix 1 — Chapter 1 Figure S1.7: Community-pair differential for combinations from three to six species simultaneously, in the four experimental treatments (control, nitrogen, fungicide, and combined), for the years (a) 2017 and (b) 2018. Asterisks above boxes indicate significant differences (p < 0.05) between the corresponding treatment and the control (Table S1.18—S1.19. 128 Figure S1.8: Community-pair differential for combinations from three to six species simultaneously, in the four experimental treatments (control, nitrogen, fungicide, and combined), calculated using the coefficients estimated for the (a) lower and (b) upper limits of the 95% confidence interval. Asterisks above boxes indicate significant differences (p < 0.05) between the corresponding treatment and the control (Table S1.20—S1.21. 129 Appendix 1 — Chapter 1 Figure S1.9: Feasible combinations (% of the total) of two to six species, in the four experimental treatments (control, nitrogen, fungicide, and combined), for the years (a) 2017 and (b) 2018. A feasible combination is one in which all species are predicted to coexist. Asterisks above boxes indicate significant differences (p < 0.05) between the corresponding treatment and the control (Table S1.23—S1.24). The total number of possible combinations (combos) of two-six species is shown at the bottom. 130 Figure S1.10: Feasible combinations (% of the total) of two to six species, in the four experimental treatments (control, nitrogen, fungicide, and combined), calculated using the coefficients estimated for the (a) lower and (b) upper limits of the 95% confidence interval. A feasible combination is one in which all species are predicted to coexist. Asterisks above boxes indicate significant differences (p < 0.05) between the corresponding treatment and the control (Table S1.25—S1.26). The total number of possible combinations (combos) of two-six species is shown at the bottom. 131