Species richness increases fitness differences, but does not affect niche differences

A key question in ecology is what limits species richness. Modern coexistence theory presents the persistence of species as a balance between niche differences and fitness differences that favor and hamper 3 coexistence, respectively. With most applications focusing on species pairs, however, we know little about if and how this balance changes with species richness. Here, we apply recently developed definitions of niche and fitness differences, based on invasion analysis, to multi-species communities. We present the first 6 mathematical proof that, for invariant average interaction strengths, the average fitness difference among species increases with richness, while the average niche difference stays constant. Extensive simulations with more complex models and analyses of empirical data confirmed these mathematical results. Com9 bined, our work suggests that, as species accumulate in ecosystems, ever-increasing fitness differences will at some point exceed constant niche differences, limiting species richness. Our results contribute to a better understanding of coexistence multi-species communities. 12


Introduction
Explaining nature's biodiversity is a key challenge for science (Hutchinson, 1957). One type of approach consists of focusing on the capacity of individual species to persist through time despite occasional pruning 15 to low density or species interactions (Turelli, 1978). Modern coexistence theory is such an approach, and predicts species persistence in a community when niche differences overcome fitness differences. Niche differences measure the strength of negative frequency dependence, i.e. whether a focal species i can recover 18 when reduced to low abundance (Adler et al., 2007;Chesson, 2000;Spaak et al., 2020). Fitness differences measure the intrinsic strength of that species (Adler et al., 2007;Hart et al., 2018;Spaak et al., 2020).
However, few applications of coexistence theory focused explicitly on explaining the persistence of 21 species in species-rich communities. Instead, most applications considered two-species communities (Chesson, 2000;Letten et al., 2017), using a variety of approaches and case studies, including annual and perennial plants (Adler et al., 2018;Godoy & Levine, 2014), phytoplankton (Gallego et al., 2019;Narwani et al., 24 2013) and bacteria (Zhao et al., 2016), evaluated under different environmental conditions (Bimler et al., 2018;Cardinaux et al., 2018;Grainger et al., 2019;Lanuza et al., 2018;Matías et al., 2018;Napier et al., 2016;Wainwright et al., 2019). To our knowledge only three empirical studies report niche and fitness 27 differences in communities composed of more than two species (hereafter multi-species communities) (Chu & Adler, 2015;Petry et al., 2018;Veresoglou et al., 2018). However, none of these three studies explained how niche and fitness differences change with species richness. In order to understand how niche and fitness 30 differences co-determine species persistence in multi-species communities, we need to understand how both variables change when adding species to a community.
Multi-species communities possess at least four complexities that are absent from two-species commu-33 nities, which may affect niche and fitness differences, and therefore how we interpret coexistence in multispecies communities. (1) First, a multi-species community can host more interaction types than species pairs, e.g. competitive or mutualistic interactions. Species richness increases the number of possible in-36 teractions and the number of combinations of these interaction types. Several metrics exist to summarize this diversity of interaction types and study their implications for community dynamics (Fontaine et al., 3 2011). (2) Second, two-species communities are always fully connected (Carpentier et al., 2021) and cor-39 relations between interspecific interactions (Barabás et al., 2016) become irrelevant since there is only a single pair of interspecific interactions. In contrast, in an n-species community there may be anywhere from n − 1 to n 2 (n − 1) connections, and the interspecific effects of species j on species i can be positively or 42 negatively correlated with the effects i has on j ( Barabás et al., 2016). May (1972) and Allesina & Tang (2012 have shown that connectance and correlation can have large effects on the local stability of multi-species communities. We therefore expect these factors to influence coexistence as well. (3) Third, 45 higher-order interactions, through which a third species changes the interaction between two species, are usually considered absent from two-species communities (but see Letten & Stouffer (2019); Levine et al.
(2017) who include higher-order interactions into two-species communities). Such higher-order interactions 48 have been found empirically, for example, in communities composed of phytoplankton, bacteria, and ciliates (Mickalide & Kuehn, 2019). In that study, bacteria coexisted with phytoplankton and ciliates, but all three functional groups did not coexist. (4) Fourth, indirect effects, whereby a third species changes the dynamics 51 of a species pair by directly interacting with one or both partners, are by definition absent from two-species communities (Walsh, 2013). We will refer to these four complexities throughout the text with (1) interaction types, (2) interaction matrix structure, (3) higher-order interactions and (4) indirect interactions.

54
Studying multi-species coexistence is challenging both theoretically and experimentally. Theoretically speaking, the methods to analyse coexistence via niche and fitness differences in a multi-species community were not available until recently (Carmel et al., 2017;Carroll et al., 2011;Chesson, 2003;Saavedra et al., 57 2017; Spaak et al., 2020). Experimentally speaking, studying coexistence of multiple species is resourcedemanding. For instance, in the simple case of linear direct interactions among species (i.e. as in Lotka-Volterra models) the number of experiments needed to parametrize the community is quadratic in species 60 richness (but see Maynard et al. (2020)). Considering higher-order interactions will consequently result in an even higher experimental load. For example, measuring higher-order interactions, sensu Letten & Stouffer (2019), would require 39 experiments in a three species community.

63
Here, we investigate the balance between niche and fitness differences along a gradient of species richness. More specifically, we ask how niche and fitness differences change as the number of species increases 4 in randomly assembled communities, and how the additional complexities (1)-(4) influence these changes. 66 We do so using four independent methods that rely on a novel, species-specific definition for niche and fitness differences (Spaak et al., 2020). First, we derive equations that quantify how niche and fitness differences respond to species richness in a community with linear interactions and simple cases of higher-69 order interactions. Second, we give an intuitive explanation of these responses based on the Mac-Arthur consumer-resource model. Third, we perform simulations with more complex models. We run these simulations as a full-factorial virtual experiment, varying direct interactions (type, correlation, connectance), 72 higher-order interactions and indirect interactions. Fourth, we compile data from the literature on empirically measured species interaction matrices and compute niche and fitness differences. We then compare our results obtained via random species assembly to a community with non-random assembly. All methods 75 support the same general conclusion: species richness does not affect niche differences, but increases fitness differences. Importantly, these conclusions are independent of the four complexities (1)-(4).

Model assumptions
To include the additional complexities of multi-species communities we use a generalized Lotka-Voltera model with n species containing second-order interactions to model the per-capita growth rates f i (N i , N −i ): where N i is the density of the focal species i, N −i is the density of all species except the focal one (vector of length n − 1), r i is its intrinsic growth rate, and α i j and β i jk are first (or linear) and second-order species interactions, respectively. We focus on basal species, e.g. primary producers, and therefore assume r i > 0.

84
A positive α i j indicates a positive interaction between species i and j (facilitation). Negative α i j indicate negative interactions (competition). Additionally, two species i and j do not necessarily have the same effect on each other as α i j can differ from α ji in both sign and magnitude. When β i jk is positive or negative, species 87 k will intensify or weaken the relationship between species i and j, respectively (second-order interaction).
The inclusion of third-order interactions did not affect any of our results. Throughout the manuscript we assume α ii = −1, which can be achieved by re-scaling of empirical data. Furthermore, we initially assume 90 a random community assembly (but see section Non-random community assembly). That is, α i j and α ik ( j = k) are sampled from any independent distributions (Appendix S1, similarly, the β i jk are all sampled from any independent distributions for k = j = i (we chose β i jk ∈ [−0.05, 0.05]).

Niche and fitness differences
There exist a total of eleven definitions for niche and fitness differences with various assumptions about the community model (Spaak et al., 2020). Five of these definitions apply to multi-species communities 96 (Carmel et al., 2017;Carroll et al., 2011;Chesson, 2003;Saavedra et al., 2017;Spaak et al., 2020). While these definitions highlight different aspects of a community model, they all capture the idea that higher niche differences facilitate coexistence, while higher fitness differences hamper coexistence. Optimally, we would 99 compute niche and fitness differences according to all five methods to understand how species richness affects coexistence.
However, the available definitions all are based on different understanding of what niche and fitness 102 differences should intuitively mean. We have shown that these different intuitions make it hard to compare across definitions. Here, we focus on one specific definition that we feel is best suited to address our specific research question (Spaak et al. (2020), rationale explained in Appendix S2). Yet, we acknowledge that 105 other interpretations of the concepts of N i and F i (i.e. different definitions of N i and F i ) will lead to different answers on how species richness affects N i and F i . To illustrate this variation, we have applied four different methods to competitive Lotka-Volterra community models (Appendix S2).

108
Importantly, the definition of Spaak et al. (2020) for niche and fitness differences capture the intuitive concepts of niche and fitness differences known from the Lokta-Volterra community model (Chesson, 1990;Chesson & Kuang, 2008). Specifically, negative niche differences imply positive frequency dependence 111 (Ke & Letten, 2018), niche differences exceeding 1 imply facilitation (Koffel et al., 2021) and intermediate niche differences imply negative frequency dependence. Stronger intrinsic strength is associated with lower 6 fitness differences (more positive F i ) and species persist, i.e. have a positive invasion growth rate, if niche differences overcome fitness differences (−F i > N i 1−N i ). Additionally, these definitions align with the niche and fitness differences for a two-species Lotka-Volterra community model Spaak et al. (2020). Spaak et al. (2020) base their definition of niche and fitness differences for a focal species i (N i and 117 F i ) on the comparison of growth rates in various conditions. If the two species (i and j) of a two-species community have completely separated niches (N i = N j = 1), the species i will grow in absence of j as it would in its presence. Which implies that intrinsic growth rate of species i ( f i (0, 0)) is equal to the invasion 120 growth rate, the per-capita growth rate of species i when it invades a community with j at equilibrium Carroll et al. (2011);Chesson (1994); Ellner et al. (2018) Conversely, if the two species have exactly the same niche (N i = N j = 0) they have equivalent effects on 123 each other. It then holds that f i (0, N * j ) = f i (c i j N * j , 0) where c i j is the conversion factor allowing to express individuals of species j as individuals of species i.
Unfortunately, c i j are defined implicitly as the solution of an equation, which often does not have a 126 closed form solution, but can be computed numerically. However, intuitively the c i j are simple. The conversion factors c i j change the frequency of species i and j, but do not affect the total density. That is, two communities with densities (0, N * j ) and (c i j N * j , 0) have the same total density, but different frequencies of 129 species i. In the first community species i has 0 frequency, while in the second it has 100% frequency, as c i j N * j is a density of species i. We refer to Spaak et al. (2020) and the appendix S3 for a mathematical definition of c i j .

132
Further intuitive insight can be gained from a Lotka-Volterra model in which the c i j = α i j α j j α ji α ii (Appendix S4). Note that the arrangement of the interaction coefficients α i j differs from both, the niche and the fitness differences. The numerator (α i j α j j ) is the competitive effect of species j, while the denominator 135 (α ji α ii ) is the competitive effect os species i. Therefore, c i j keeps total density constant by scaling the effects of the two species.
Niche differences are then defined via interpolation between these two extreme cases of N i = 0 and 138 7 N i = 1: This equation of niche differences represents the relative frequency dependence of species i. The numerator assesses frequency dependence, as it compares the growth rate of species i when it has 0% frequency ( f i (0, N * j )) to when it has 100% frequency ( f i (c i j N * j , 0)), notably at the same converted density. The denom-

144
inator assesses density dependence, as it compares the growth rate of species i at 0 density ( f i (0, 0)) to when it has non-zero density ( f i (c i j N * j , 0)). Importantly, in both growth rates species i has 100% frequency.
Similarly, Spaak et al. (2020) define fitness differences as the scaled growth rate in the absence of niche Zero fitness differences imply that the species have equal intrinsic strength, as N * i = c i j N * j . That is, species will have the same equilibrium density after conversion with c i j . Competitive subordinate species then have negative fitness differences, as N * i < c i j N * j , conversely competitive dominant species have positive fitness 153 differences. Note that F i = 0 means no fitness differences and more negative F i means stronger fitness differences. Importantly, these definitions of niche and fitness differences agree with the original definitions on the two species Lotka-Volterra models (Chesson, 2013;Chesson & Kuang, 2008), they are therefore not 156 newly defined niche and fitness differences, but rather a generalization thereof (Appendix S4, Spaak et al. (2020)).
Additionally, this definition, illustrated so far for a two-species community, naturally generalizes to 159 multi-species communities. The three growth rates remain conceptually the same. The intrinsic growth rate

Analyses and Simulations
We first examined analytically how N i and F i change with species richness. We found a generic solution for first-order interactions and for a simplified case of higher-order interactions.
Second, we designed a full-factorial virtual experiment in which we numerically computed N i and F i for a wider range of different simulated communities (Table 1). For these we solve numerically for equilibrium densities and invasion growth rates using the 'fsolve' function from the scipy package in Python.

174
Communities with higher-order interactions can have multiple equilibria (Aladwani & Saavedra, 2019). To compare communities with and without higher-order interactions we computed N i and F i only for the equilibrium which was closest to the equilibrium without higher-order interactions (Appendix S1). We did so 177 because the chosen higher-order interactions are small, i.e. β i jk ∈ [−0.05, 0.05], such that the higher-order interactions should be seen as a perturbation of the linear interaction. The factors used in our simulations were (i) first-order interaction type (competitive, facilitative or both, i.e. α i j < 0, > 0 or unrestricted); (ii) 180 first-order interaction strength (strong or weak); (iii) connectance of the interspecific interaction matrix(c ∈ 1, 4 5 , 2 3 ); (iv) correlation between the interspecific interaction (ρ(α i j , α ji ) = ρ i j (β i jk , β jik ) ∈ {−1, 0, 1}); (v) inclusion of indirect effects (present or absent); (vi) second-order interaction type (β i jk < 0, > 0, un-183 restricted, or absent). To exclude indirect effects we set equilibrium densities of resident species to their monoculture equilibrium density. In this way, we cancel out interactions among residents that will change the residents' densities. The intrinsic growth rate r i does not affect N i and F i , therefore we set it to 1 in all 186 simulations. However, the r i may affect the community dynamics from stable coexistence to chaotic fluctuations and competitive exclusion (Song et al., 2020). In this case, invasion growth rates do not correctly predict coexistence and should not be analyzed with N i and F i (see limitations).

189
This design leads to a total of 3 · 2 · 3 · 3 · 2 · 4 = 432 parameter settings. We ran 1000 repetitions for each of the five species richness levels (2 ≤ n ≤ 6), leading to a total of 432 · 5 · 1000 = 2 160 000 simulations. We parametrized the first-order interactions using distributions of empirically obtained first-order 192 9 interactions (Appendix S1). We sampled "strong" first-order interactions between the Q 1 − 1.5(Q 3 − Q 1 ) and Q 3 + 1.5(Q 3 − Q 1 ) of those distributions, where Q 1 and Q 3 are the first and third quartile, to remove outliers (Appendix S1). Similarly, we sampled "weak" first-order interactions between the Q 1 and Q 3 of the 195 empirical distributions of interaction strength. We fitted linear regressions to measure the effect of species richness on N i and F i . With this approach we were able to investigate the effects of all complexities (1)-(4).
Literature data 198 We found three review papers that collected multi-species Lotka-Volterra interaction matrices (Adler et al., 2018;Fort & Segura, 2018;Keddy & Shipley, 1989), representing a total of 33 interaction matrices, ranging from 3 to 9 species, and containing 29 plant, 2 phytoplankton, 1 zooplankton and 1 ciliate communities. We 201 normalized all these data such that α ii = −1. The interaction matrices were obtained through pairwise experiments, measuring the interspecific effect of one species on the other. For each multi-species community we constructed all possible sub-communities with at least two species, leading to a total of 2544 commu-204 nities that varied in species richness from 2 to 9. We excluded all communities in which not all interaction strengths were available, e.g. because of a "NA" entry in the sampled sub-community, leading to 2296 communities. For 1376 communities we could not compute N i and F i . That is because, like any method 207 seeking to quantify frequency dependence, our approach is based on invasion analysis: the capacity of an invader to grow with the other species at their non-zero equilibrium. N i and F i are thus only computable for communities where all species in each subcommunity (the community without the invading species) coexist 210 stably. We therefore computed N i and F i for the remaining 920 communities (Appendix S5).
For each interaction matrix obtained from the literature we computed N i and F i using equation 2 and 3.
For each of the 33 interaction matrices, we regressed N i against species richness of the sub-communities. 213 We used a Theil-Sen estimator for the slope, which is more robust to outliers than linear regression based on least squares (Sen, 1968). We fitted (using least squares) a saturating function F i = n−2 (n−2)+H for the fitness differences, where n is the species richness and H is the free parameter. This saturating response was chosen 216 for F i , because our analytical results suggested such a response.

Analytical solutions 219
For the linear Lotka-Volterra model without higher-order interactions (i.e. β i jk = 0), we can explicitly compute the fitness F i and niche differences N i of species i in the multispecies community (Appendix S4): with F (2) i j and N (2) i j the fitness and niche difference of species i in a two-species community consisting of 222 species i and j. N −i, * j is the equilibrium density of species j in the absence of species i, N * j is the equilibrium density of species j in monoculture and c i j is the conversion factor from species j to species i (Methods).
The sums are taken across all species j with which species i interacts directly, i.e. α i j = 0.

225
Equation 4 and 5 illustrate two main results. First, 1 − F i is the weighted sum of the two-species fitness differences 1 − F i j , across all species pairs with which i interacts. The weights are the relative yields, as known from biodiversity ecosystem functioning research (Fox, 2005;Hector & Loreau, 2001). The effect 228 of species richness on fitness differences will therefore be similar to the effect of species richness on the sum of the relative yield, known as the relative yield total (∑ j =i . The relative yield total is known to increase with species richness in many communities (Carroll et al., 2011;Grace et al., 2016;Loreau, 231 2004). Hence, 1 − F i will, on average, increase with species richness. Second, N i is the weighted average of the two-species niche differences N i j . Hence, species richness will, on average, not affect N i . Since we did not make assumptions about the distributions of α i j , these results are independent of the details of 234 interspecific interactions, i.e. the results apply regardless of complexities (1) and (2).
Equation 4 and 5 can be approximated by assuming constant interspecific interaction strength α (Appendix 4). This yields N i ≈ 1 − α and F i ≈ 1 − n−1 1−(n−2)ᾱ , from which it is clear that N i is independent of because the sum of the weights , the relative yield total, will saturate as well in the Lotka-Volterra model (Loreau, 2004;Spaak et al., 2017).

240
Link to resource competition The fact that 1 − F i is a weighted sum, while N i is a weighted average makes intuitive sense when realising that the interaction coefficients α i j can, under certain conditions, be related to resource utilisation (Chesson, 243 1990;MacArthur, 1970). Consider a focal species (yellow, Fig. 1) utilizing a range of resources. We ask how adding competitors changes this species N i and F i (Fig. 1 A-E). We assume that the species only differ in resource utilisation rates, not in other parameters such as mortality.

246
In a two species community, the species with the higher total utilisation rate will have a competitive advantage and consequently the higher fitness difference. The total resource utilisation rate, denoted A i , is the area under the curve A i . One could therefore intuit that the fitness difference is linked to the ratio 249 of total resource utilisation rates, i.e. F 1 ≈ 1 − A 2 A 1 . Fitness differences therefore increase with species richness, as each competitor increases the total resource utilisation rates of all competitors combined ( Fig.   1 F), i.e. F 1 ≈ 1 − ∑ A j A 1 . It turns out that this intuition is almost correct; we only have to add weights to 252 the sum according to the densities of the species at equilibrium (compare this equation to eq. 4). F 1 thus increases (becomes more negative), as species richness increases.
Intuitively, 1−N i is the proportion of the shared resources between the focal species and its competitors, 255 that is the amount of the shared resources scaled by the total consumption of the species. In a two species community, we therefore intuit N 1 ≈ 1 − A 1 ∩A 2 A 1 · A 2 . Increasing species richness will increase the sum of the shared resources, as the focal species will share resources with each competitor (Fig. 1 G), but also 258 the sum of the total consumption of the species (Fig. 1 F). We therefore expect N 1 ≈ 1 − A 1 ∩∑ j A j A 1 ·∑ j A j to be independent of species richness (Fig 1 H). Again, this intuition is correct up to the inclusion of the weights according to the species equilibrium densities.

Full-factorial simulations
The simulations with strong first-order interactions only partially seem to confirm the predictions made by theory (Fig. 2). That is, N i is not invariant but approaches 1 as species richness increases (Fig. 2 A), which 264 seems to contradict the theoretical results. Yet, species richness does not directly affect N i , but rather affects the average interaction strengthᾱ, which in turn affects N i , which can be approximated by N i ≈ 1 −ᾱ (Appendix S4). That is because, by design, any method based on invasion growth rates (such as those to 267 compute N i and F i ) can only be applied to communities in which invasion analysis is possible. Hence, too strong negative interactions prevent the invasion into highly-diverse communities, and will often impede feasible n − 1 subcommunities to begin with (Kokkoris et al., 2002). Hence, species richness selects against 270 communities with overly strong negative interactions, which leads to, on average, weaker interactions at higher species richness (Appendix Fig.S3). Similarly, species densities in communities with strong positive interactions will tend to grow to infinity, and more so in species-rich communities, because interspecific 273 facilitation is stronger than intraspecific limitation (self-regulation). Again, species richness selects against strong positive interactions, weakening the average interaction strength (Appendix Fig.S3). This selection of weak (negative and positive) interspecific interactions causes N i to approach 1. F i increases with species 276 richness for all parameter settings, as predicted by the theory (Fig 2 B).
The simulation results based on weak interaction strengths allow assessing the direct effect of species richness on N i and F i without the confounding effect of species richness on interspecific interaction strength 279 α i j . In these simulations, the effect of species richness on interspecific interactions was much weaker (Appendix Fig.S3). These simulations confirmed our theoretical results; N i was on average unaffected by species richness (Fig. 3 A) and F i increased with species richness (Fig. 3B). We illustrate how N i and 282 F i values jointly varied with species richness, using weak interaction strength: no higher-order interactions (β i jk = 0), no correlation between the α i j , and maximum connectance (Fig. 3 C). Again, these results hold independently of the complexities (Appendix S1).

285
As expected, using other definitions of N i and F i to the simplest case (competitive Lotka-Volterra communities) leads to a variety of responses of N i and F i to species richness. Some of these definitions 13 respond in the same way to richness (e.g. Chesson (2003) for F i ), while others respond in an differently 288 (e.g. (Carroll et al., 2011) for N i ). This comparison shows that it is crucial to use the appropriate definition of N i and F i for the question at hand (Appendix S2).

291
The results for the empirical communities reflect those obtained for the simulated communities. The absolute values of the slopes of the linear regressions of N i were small (< 0.05) for all but 6 datasets. The slope for the overall regression of N i against species richness (Fig. 4A, black line) was small (-0.028). F i increased 294 with richness in all but one dataset. Overall, we conclude that the response of N i and F i to richness for empirical communities did not qualitatively differ from that of randomly generated communities.
The empirical data also revealed cases in which coexistence is possible even though some of the species 297 have negative N i . This is possible as long as F i is sufficiently positive such that F i ≥ −N i 1−N i . A total of 95 (4.1%) communities were found with species persisting despite having negative N i .
Non-random community assembly 300 So far we have focused on random community assembly, where the location and width of the resource utilization A i of species i were chosen randomly ( Fig. 1 A-E). We here expand our results towards two different possibilities of non-random community assembly. First, given a species pool, we rearranged the 303 order of species arrival. By choosing a non-random community assembly, species richness can increase or decrease niche differences, but will always increase fitness differences (Appendix S6). Additionally, averaged over all possible community assemblages, species richness again does not affect niche differences.

306
Alternatively, one might ask how species richness affects niche and fitness differences, if species optimize their resource utilisation according to the prevailing species richness (Barabás et al. (2016), Fig. 5 A-E, Appendix S7). Importantly, in this case, the two-species community is not a subcommunity of the three-309 species community. This links to the traditional question of species packing, which is how close species can be packed in a given environment and still coexist (Barabás et al., 2014;MacArthur, 1970). In this scenario, species richness still increases fitness differences (Fig. 5 F). However, because of the limited niche space, species in species rich communities are located closer to each other and therefore have a stronger interaction strength on average (Fig. 5 G). Niche differences therefore decrease with species richness (Fig. 5 H).

315
It is well-established that the likelihood for stable coexistence drops with species richness Goh & Jennings, 1977;May, 1972;Serván et al., 2018). Here, we explain this well known result in the context of modern coexistence theory by examining how niche and fitness differences (N i and F i ) 318 change with species richness. We found that species richness has no direct effect on N i but directly increases F i . However, species richness may have an indirect effect on N i (Figure 2 and 5). This conclusion is based on four independent approaches: mathematical computation, biological intuition, numerical simulations, 321 and analysis of experimental data. Overall, the influence of species richness on N i and F i is robust to inclusion or omission of the complexities (1)-(4), and all their combinations. The fitness differences of a species increases with species richness, as fitness differences measure the fitness of a species compared to the com-

324
bined fitness of all other species. In multi-species communities, most species will therefore have negative fitness differences, as one species will rarely have higher fitness than all other species combined. Moreover, N i in multi-species communities is approximately 1 −ᾱ, whereᾱ is the average inter-specific interac-327 tion strength. This new result facilitates considerably the computation of N i in multi-species communities.
Taken together, our findings shed new light on the causes of coexistence in multi-species communities.
The niche differences of a species measure the proportion of limiting factors, e.g. resources, that are 330 limiting to other species as well. Increasing species richness increases the amount of limiting factors shared with other species, but also the amount of limiting factors that are not shared with other species. The proportion of shared limiting factors is therefore unaffected, on average. Species-rich communities are therefore less likely to coexist (all else being equal), as fitness differences become too strong to be overcome by niche differences.

336
The available experimental data only represent fully connected communities, with no correlation among interactions (complexity (2)) and, most notably, did not contain cases of higher-order interactions (complexity (3)). We do therefore not know whether the parameter values used to describe these higher-order interac-339 tions in our simulations (and therefore the simulation results) are realistic. The available experimental data were biased towards competitive communities of terrestrial plants with relatively low species richness. Our simulations suggest that our conclusions hold for other networks as well. Computing N i and F i on a larger 342 collection of natural communities would help refine our understanding of this process. However, obtaining the full interaction matrix for species-rich communities is challenging. To obtain interaction matrices, various approaches exist. For example, one can use the frequency of interaction between species (e.g. number

345
of visits of a pollinator on a plant) as a proxy for interaction strength. The robustness of this approach, however, still needs to be tested (García-Callejas et al., 2018). Other methods rely, for example, biomass (Moore et al., 1996;Zhao et al., 2019), mass ratio (Emmerson & Raffaelli, 2004) or production and consumption 348 rates of species (Christensen V. & D., 1992;Jacquet et al., 2016). These different methods rely on different assumptions and may therefore influence the resulting matrix estimate (Carrara et al., 2015).
As mentioned before, the computation of niche and fitness differences is based on invasion growth rates, 351 which may not always be defined (Barabás et al., 2018;Spaak et al., 2020) or may not predict coexistence (Pande et al., 2019). This makes it especially challenging to assess the importance of indirect species interactions, and more specifically, intransitivity. For example, a community coexisting solely through intransitive 354 coexistence could not be analyzed with our method as invasion analysis will not be possible.
Given these limitations, one can ask to what extent addressing them would change our conclusions.
In communities where species richness increases total abundance, which is the case for various com-357 munities (Grace et al., 2016;Loreau, 2004;Turnbull et al., 2013), we expect the no-niche growth rate , 0) to become more negative, as ∑ j c i j N −i, * j increases (eq. 2, 3). Consequently, we expect species richness to increase fitness differences, i.e. make F i more negative, as it is linear in the no-niche 360 growth rate. Conversely, in communities where species richness decreases total abundance we expect the opposite, that is: fitness differences might decrease with species richness. It is less clear how species richness will affect niche differences in models not explored in the current paper , e.g. with different per-capita 363 growth rates functions f i , or with a community with age structure (Chu & Adler, 2015). Niche differences depend on the invasion growth rate and the no-niche growth rates, which both depend on the species richness and total abundance. When species richness has a stronger negative effect on the no-niche growth rate than 366 on the invasion growth rate, then niche differences will increase with species richness. Conversely, if the invasion growth rate decreases more than the no-niche growth rate, niche differences will decrease.
Additionally, we have mostly assumed that species richness is independent of other factors such as the 369 average interaction strength, connectance or correlation (Cohen & Briand, 1984;Kokkoris et al., 2002). It is well-established that the number of links scales nonlinearly with species richness, causing connectance to decrease with species richness (Carpentier et al., 2021;Cohen & Briand, 1984;Martinez, 1992). However, 372 this will not affect our findings, as the sums in the equations 4 and 5 only run over species with which the focal species i interacts. We ran additional simulations to illustrate this point (Appendix S8). As predicted by theory, niche differences are unaffected and fitness differences increase with species richness. However, 375 fitness differences increase more slowly than in the other investigated cases, as the number of links (and therefore the number of summands in equation 3) increases more slowly.
The main determinant of niche differences is the average interaction strength, which might decrease 378 with species richness due to coexistence requirements ( Fig. 2A) or might increase due to increasing overlap of resource requirements (Fig. 5 G). How species richness covaries with average interaction strength will affect how species richness covaries with niche differences, it is therefore very likely that species richness 381 and niche differences will covary in natural communities.

New insights
Our results yield two new insights, other than the main result on how N i and F i varies with species richness.

384
A first insight is that negative niche differences do not necessarily preclude coexistence. Negative niche differences have been attributed to priority effects and therefore viewed as precluding coexistence (Fukami et al., 2016;Ke & Letten, 2018). Our framework confirms this finding for the case of competitive two-species communities (Spaak et al., 2020). However, species in multi-species communities will not all have the same niche differences (example community in Fig 4). This implies that species a, with negative niche differences and low fitness differences, can coexist with species b and c that have positive niche differences 390 and negative fitness differences. Consequently, multiple species can have negative niche differences in a multi-species communities and still persist. In our empirical data set, we found six three-species communities in which all but one species had negative niche differences.

393
A second insight is that one can infer N i and F i in multi-species communities from N i and F i measured in pairwise interaction experiments. If one measures N i and F i for each two-species sub-community of an n species community, which is typically done ( (2) i j (n−1) . With one additional multi-species experiment to estimate the relative yield RY i we obtain an estimation of i j ) · RY j as well. This indicates that two-species experiments are sufficient to predict N i and F i in multi-species communities. To 399 validate this finding one would optimally conduct all the necessary n multi-species experiments to measure the invasion growth rates, and fit a community model to the multi-species communities. Given these invasion growth rates and the community model one can compute exact niche and fitness differences and compare 402 them to their estimates.
Finally, one of the key questions in community ecology is whether niche differences are strong enough to overcome fitness differences and allow coexistence. Often, niche differences are found to be not only 405 sufficiently strong, but much stronger than strictly needed (Chu & Adler, 2015;Levine & HilleRisLambers, 2009). The present results offer a potential explanation for this observation. That is, niche differences not only need to be sufficiently strong to overcome fitness differences of one or few competitors, as typically 408 considered in empirical studies, but sufficiently strong to overcome fitness differences of the entire resident community, as niche differences are independent of species richness. Our results therefore allow asking the more general question of how many species one can pack in a community, given its niche difference.   Figure 1: A-E: Resource utilisation of the yellow focal species and its competitors. F: Increasing the species richness will increase the total utilisation of the resident species ∑ j A j . Similarly, we expect to increase with species richness, as it scales with the ratio of total resource consumption.
G: The amount of shared resources (hatched region from panels A-E) increases with species richness. H: As both the amount of shared resources increase (panel G) and the total utilisation (panel F) increase, we expect the ratio to be independent of species richness. Similarly, we expect N i ≈ 1 − unconstrained interspecific interactions are shown in blue). However, this is not a direct effect of species richness on niche differences, but rather an indirect effect through a decreasing average interaction strength with increasing species richness (Fig. S2). Each line represents a linear regression of niche differences as a function of species richness for one factorial setting of the full-factorial experiment (Table 1). B: Species richness, however, makes fitness differences more negative (i.e. larger). Note the differences in y-scale between panel A and B. C: Distribution of N i and F i for simulated theoretical communities that are fully connected, and exhibit first-order interactions without correlations, i.e. similar to the experimental communities (Fig. 4). Each dot represents N i and F i of one species in a community composed of 2-6 species (see colour legend). The black line indicates the persistence line, species below this line are assumed to persist in the community. Note the inverted y-axis. i.e. species having a positive net effect on another, and therefore N i > 1 is common in the datasets we found. We highlight one specific three-species community (grey line) where all species coexist, even though species a has N i < 0, a property associated with priority effects and therefore exclusion. Axis from C are truncated to show ∼95% of all data points. 29  Barabás et al. (2016). F: Increasing species richness increases fitness differences, as predicted by theory G: As the niche width is limited, species are located closer to each other in species rich communities and therefore have a stronger interaction strength on average. Each dot represents one value of the interaction matrix α i j . H: Consequently, niche differences decrease with species richness, contrary to theory. The example deviates from the theoretical predictions, as one of the key hypothesis of random assembly is not meat.