The effect of natural enemies on the coexistence of competing species - an empirical test using Bayesian modern coexistence theory

The role of natural enemies in promoting coexistence of competing species has generated substantial debate. Modern coexistence theory provides a detailed framework to investigate this topic, but there have been remarkably few empirical applications to non-plant systems. Trade-offs between defence against natural enemies and inherent fecundity provide a potential mechanism promoting coexistence between competing species. To test this, we parameterised models of competition between six Drosophila species in experimental mesocosms, with and without a generalist pupal parasitoid. We found no evidence of a general fecundity-susceptibility trade off, and idiosyncratic impacts of parasitism on pairwise coexistence. Methodologically, our novel Bayesian approach highlights issues with the separability of model parameters within modern coexistence theory and shows how using the full posterior parameter distribution improves inferences. Our results emphasise the importance of contextualising specific trade-offs and the value of modern coexistence theory in multi-trophic contexts.


Introduction
Species compete for limited resources and (in almost all cases) are consumed by other species. Understanding how the impacts of competitors and natural enemies interact to influence species coexistence is a long-standing challenge. Natural enemies are regularly cited as a driver of coexistence, for example by suppressing otherwise-dominant species (Paine 1966) but their impact can be highly variable (Chase et al. 2002).
Analyses on the effect of natural enemies often focus on 'reduction in competition' (e.g. Gurevitch et al. 2000), as measured by negative effects between focal species, but consumers can affect competition and coexistence in convoluted ways (Chase et al. 2002;Pringle et al. 2019). Modern coexistence theory (MCT) focuses on the property of mutual invasibility, and allows the precise framing of questions regarding pairwise species coexistence as a balance between niche differences promoting stabilisation and fitness differences fostering exclusion of one species by another (Chesson 2000;Letten et al. 2017). Through the lens of coexistence theory, interactions mediated through competition for resources and through natural enemies are conceptually equivalent (Chesson & Kuang 2008;Chesson 2018) and a body of theory has developed to integrate consumer effects within the coexistence framework (Chesson & Kuang 2008;McPeek 2019). This theoretical work has emphasised that impacts of any process on coexistence can only be determined with information on both intrinsic growth rates and competitive terms. However, a corresponding body of direct empirical data using the framework to investigate the impact of consumers on coexistence is yet to develop .
Natural enemies are theorised to enhance coexistence through a number of mechanisms. A frequently-invoked mechanisms is that they reduce fitness differences, mediated by a trade-off between resource acquisition and susceptibility to attack (Paine 1966;Kraaijeveld & Godfray 1997). However, such a trade-off is in itself insufficient to enable long-term coexistence, since some kind of stabilising niche-difference is required to mitigate any fitness difference (Chesson 2000). Even with zero fitness difference, stochastic drift would preclude coexistence in the absence of stabilising mechanisms (Adler et al. 2007).
The great majority of experimental work within a formal MCT framework has been conducted on plant or microbial systems Grainger et al. 2019), but experimental studies of insect communities have considerable potential to further our understanding of the role of natural enemies in species coexistence. While plants are particularly amenable for competition studies, it is challenging to explore how herbivores affect coexistence because, with the exception of seed predation (Petry et al. 2018), it is complex to quantify the effect of damage by natural enemies directly in terms of fitness. By contrast, since parasitoids kill their insect hosts as they develop in or on them, each attack is directly interpretable in terms of host fitness (Hassell 2000). Although ecologists have been competing insect species (especially Drosophila) against each other for decades (Pearl 1932;Ayala 1969;Gilpin et al. 1986;Worthen 1989;Davis et al. 1998), very few studies have explicitly tested for mutual-invasibility in insect systems (Siepielski et al. 2018). In particular, few studies have directly assessed coexistence per se, rather than particular conditions necessary for coexistence (Godwin et al. 2020;Spaak & de Laender 2020). environmental changes will create novel communities (Alexander et al. 2016). However, parameter fitting can have high data demands, and key ecological model parameters are often not fully separable (Bolker 2008). For example, in a simple model of competition similar predictions could be made by lowering the growth rate term and increasing the competition term. With a well-designed experiment and sufficient data, uncertainty can be minimised. Nonetheless, in any realistically-sized competition experiment, separability of parameters can be a significant but largely underacknowledged issue. This problem is particularly pertinent when working within the framework of modern coexistence theory, where conclusions are drawn from comparisons between ratios composed of multiple parameters from the underlying model (Barabás et al. 2018;Song et al. 2019). Failing to address these issues could dramatically affect the confidence with which conclusions are drawn from experimental data.
Here we present one of the first analyses directly quantifying the impact of consumers on pairwise species coexistence, focusing on a community of six rainforest Drosophila species. Our experiment was designed specifically to capture the key parameters of our focal species and their competitive interactions to test the hypothesis that trade-offs between intrinsic growth rate and susceptibility contribute to their coexistence. We use a Bayesian approach to propagate uncertainty in model parameters through to niche and fitness differences, a method that should be widely applicable for understanding species coexistence processes in a range of systems. We show that parasitoids strongly influence the fitness differences between species and have inconsistent impacts on pairwise coexistence that would not be identifiable solely from fecundity-susceptibility relationships.
Drosophila species have long been used in experimental ecology since they are straightforward to culture under laboratory conditions, allowing high levels of replication and standardisation. The set of Drosophila species used here represents a significant fraction of the diversity of Drosophila species co-occurring at rainforest sites in Queensland, Australia (Jeffs et al. 2020). The six species span a range of body sizes (Table S1.1), have similar generation times, and can be distinguished visually (at least in the case of males). Drosophila are generalist feeders (Valadão et al. 2019) and in natural populations they experience high levels of parasitism, including from a suite of highly generalist natural enemies (Jeffs et al. 2020).

Experimental design
We conducted a large, single-generation pairwise competition assay, in a response-surface design following the guidelines of Hart et al. (2018), under the presence or absence of pupal parasitoids. We used 25mm-diameter standard Drosophila vials (VWR International) containing 4 ml of standard fly medium. Equal numbers of male and female Drosophila were introduced to each vial, so we refer to the number of founder pairs. Each two-species combination (15 combinations in total) was represented by (3,3), (6,3), (9,3), (12,3), (3,6), (3,9) and (3, 12) sets of pairs, each replicated three times in each parasitism treatment. We also included monocultures of each species, with 12x 1 pair, 6x 3 pairs, 6x 9 pairs, and 6x 15 pairs within each treatment, making a total of 990 experimental vials (360 single-species, 630 two-species). The entire experiment was split over three blocks, each with identical sets of combinations, but staggered by two days. The experiment was conducted in a single large incubator maintained at 26°C, with the location of trays within the incubator regularly rotated. Further details are given in SI 1.

Initiation of treatments
Founding adults were offspring of adults drawn from mass-bred line cages and were raised in vials at moderate density (approximately 100 adults per vial) at a constant 23°C. These adults emerged 8-10 days before the start of the experiment, and so could be assumed to be sexually mature and mated. The day preceding the start of the experiment, the founding flies were lightly anesthetised under CO2 to be separated by sex. The males were added to the experimental vials immediately. Females were maintained at moderate density (~50 adults per vial) in holding vials with fly medium and were transferred to experimental vials the following day, minimising the time interval between the addition of each species in multi-species assays. They were able to lay eggs and re-mate for 24 hours before being removed. Fig. 1 provides an overview of the experimental protocol.

Parasitism
At the onset of pupation, 5 days after the founder flies were removed, the fly vials were moved into Perspex parasitism boxes (35x20x14cm) containing up to 88 vials (SI 1). Vial positions within each box were randomised. The vial plugs were removed, and 50 female and 50 male mature and mated parasitoid wasps were added. The wasps were free to move between vials within each box. While it is possible that late-pupating larvae could cross from one vial to another during this period, we only observed a handful of cases where 1-2 individuals of the 'wrong' species emerged from a vial, and very few pupae (<5 per box) were observed to develop on the outside of the vials. This suggests that movement was negligible. After 30 hours, before the emergence of the fastest-developing Drosophila species, all wasps were removed with an aspirator, vial plugs were replaced, and the vials were returned to their original trays.

Counting
We removed adult flies emerging from all vials every three days, minimising the risk that emerged adults started a subsequent generation. Flies in single-species vials were sexed and counted at each removal time, while flies in multi-species vials were sexed and identified at the end of the experiment (Fig. 1). We observed a highly consistent 1:1 sex ratio across the intra-specific trials and in those interspecific trials where females of different species could be reliably distinguished (SI 2 Fig. S2.1). In the five of the combinations (BIR-PAN, BIR-PSA, BIR-SIM, PSA-PAN and PAL-SUL) where it was not possible to identify female flies to species we counted the males and inferred the total number of flies of each species assuming a 1:1 sex ratio.

Determining growth rates
The larval competition of Drosophila often occurs within individual fruits that represent a discrete and temporary resource suitable for only a single generation of flies (Atkinson & Shorrocks 1981;Shorrocks 1991). In this way, the population dynamics of Drosophila resemble those for annual plants that have been the focus of the majority of recent coexistence work (Levine & HilleRisLambers 2009;Pérez-Ramos et al. 2019). Our starting point is a Beverton-Holt model without any carry-over between generations: where , is the size of the population of species in generation t, is a species-specific intrinsic growth rate term and are competition terms representing the effect of species on .
The mean number of adults emerging per founder was consistently lower in vials with only a solitary founding pair compared to those with three founders (SI 2, Fig. S2.2). Drosophila use cues from other females when determining oviposition sites (Wertheim et al. 2002;Duménil et al. 2016), which could be causing this Allee effect. It is not possible to directly incorporate such positive densitydependence effects within the standard MCT framework without additional assumptions (Barabás et al. 2018;Schreiber et al. 2019). Since solitary founders are quite possibly a rather artificial scenario, we therefore excluded data from all singleton pair trials from the model fitting.

Variable generation time
One key difference between strictly annual plants and our Drosophila system is that the time from egg-laying to emergence can be considerably extended in high-competition environments, reducing the effective rate of reproduction. Additionally, our test species differ in baseline development time; in particular D. birchii has a longer development time. To account for this, we calculated an effective daily growth rate as our core 'fitness' measure from the raw number of emerged adults over a period of time. We use a minimally simple growth model, 0 × λ = , to identify a daily growth rate λ: To avoid observed counts of zero overly influencing the dataset, we added 0.1 to all values of . This approach maintains the intuition that observing no emergences from a larger number of founders would imply a lower estimated growth rate than the same observation from fewer founders. As an alternative approach, we also fitted models directly to the observed counts. This necessarily ignores development time effects, but allows direct representation of estimation errors at low counts via a model with a negative-binomial distribution error structure. This alternative analysis, which shows broadly similar results, is presented in SI 5.
Logistical constraints meant we could only observe the mean total generation time directly in the monoculture trials. We used these to estimate the emergence time in the interspecific trials with a linear regression model selected using AIC, using as predictors the number of founder flies, the species, the treatment, and all 2-way interaction terms. This model had an adjusted 2 of 0.785 for the monoculture data. The emergence of larger fly species was delayed for longer by high densities than the emergence of smaller fly species (SI 2).

Inferring model parameters
The models were fit within a Bayesian framework and posterior distributions of the parameters were obtained from Hamiltonian Monte-Carlo sampling using STAN (Carpenter et al. 2017). Since we counted both species, each two-species experiment contributed two data points. Although this contributes a small component of non-independence, the gain in total sample size is significant.
Since we excluded the single-founder vials, our effective total 'n' across both treatments was 1476, to estimate a maximum of 84 free parameters in our most complex model.
For each vial, the inferred daily growth rate λ for a particular species in the possible presence of species was modelled as: and fit assuming a Gaussian error distribution. A separate error term was fit for each parasitism treatment. The competition α terms were fit with an underlying Gaussian prior of mean 0 and σ 1 and constrained to be positive, implying only competitive interspecific interactions, with no facilitation possible. Growth rate terms were constrained to be positive, and fit with a weak Gaussian prior with mean 1 and σ 1. Tests fitting α values with different priors did not give notably different results (SI 3).
We fit and compared a sequence of candidate models of increasing simplicity. In the first, separate and α terms were fit for each treatment. In the second, a single set of α terms, but separate terms were fit across both treatments, while in a third model a single set of parameters were fit across both treatments (i.e. parasitism has no effect). The three models were compared based on expected log pointwise predictive density (ELPD, a measure of out-of-sample predictive capacity), computed using the loo R package (Vehtari et al. 2017). This approach allows an assessment of whether there is sufficient evidence in the data to separate the parameters assigned to the treatment groups in a comparable manner to other information criteria such as AIC (Burnham & Anderson 2002). All STAN model code and data to fit the models are available online. Posterior predictive checks and diagnostics are presented in SI 4.

Identifying Niche and Fitness differences
Following Chesson (2000) and Godoy & Levine (2014), we define ρ as the niche overlap between two species, and 1 − ρ as the niche difference. For species and this is calculated as: In certain combinations of α values, niche overlap ρ calculated in this manner can be greater than one, implying a 'negative niche difference'. This does not have a clear ecological meaning, and some previous authors have elected to define such cases as 0. However, since these negative 1 − values can be ecologically interpreted as potentially identifying priority effect situations (Ke & Letten 2018), we present the values as calculated. Whether negative or zero, in either case coexistence is not possible in the MCT framework.
Within the framework, fitness differences are calculated as a ratio: Coexistence is determined to be possible if: For each pair of species, we calculated the niche and fitness differences across 3000 draws from the posterior distributions of parameter estimates. We partitioned the posterior to infer the support given by the data to each of the four outcomes of competition: coexistence, victory by either species or a situation of priority effects (see Fig. 3a). Note that the dynamics described by our models are fundamentally deterministic -these values express uncertainty in outcome given the data, rather than necessarily stochastic outcomes.

Results
Based on model comparison, there was strong support for the expectation that parasitism reduces the growth rate of each species, but little support for changes in competition coefficients between the treatments. Our second model (with different values between treatments but common α) performed best (ΔELPD vs Model 1 = -10.9, vs Model 3= -51.37, Table S4.1). The parasitism treatment reduced population growth rates, but we did not identify a susceptibility -fecundity trade-off across our six species (Fig. 2).

Figure 2
No clear susceptibility -fecundity trade-off across the six Drosophila species (Spearman's rank correlation test on means: p = 0.497). Mean value of each species is plotted alongside a posterior distribution point cloud and fitted ellipse (95%, fitted assuming a multivariate-t distribution) to show distribution of uncertainty for each species. We present susceptibility in terms of change in r-1, since this is the quantity used in the calculation of fitness differences.
The impact of the parasitoids on coexistence was inconsistent across species pairs (Fig. 3c). In three species-pairs the fitness difference was so large there was essentially no support for coexistence under either treatment. Of the remaining twelve pairs, five had higher support for coexistence under the parasitism treatment, while seven had less support (Table 1). In only one pair (D. sulfurigaster and D. pandora) did the parasitism treatment cause the median estimate of fitness difference to cross the boundary delineating coexistence from competitive exclusion. However, this was associated with only a small decrease in the posterior support for coexistence (Table 1). Overall, intraspecific effects were comparable in strength to interspecific effects (SI 4, Fig. S4.2) implying a high niche overlap within many of the species pairs. This had the consequence that 'negative niche differences' (where coexistence is not possible) were identified in over half (59%) of posterior draws.
The fitted posterior distributions showed strong positive correlations between the and α values (mean Pearson correlation = 0.76), resulting in non-independence of estimates of the fitness and niche differences and non-circular posterior distributions (Fig. 3). The positive correlation between α and can effectively make the coexistence boundary line somewhat 'repelling'. In many cases, where the posterior was centred outside the region of coexistence, the support inclines away (e.g. D. sulfurigaster and D. birchii), while if the support is centred inside, the support tends to stay inside the zone of coexistence (e.g. D. simulans and D. pandora). Figure 3. The likelihood of coexistence of Drosophila species pairs responds idiosyncratically to parasitism. a) Guide to interpretation of coexistence plane diagrams. The region of coexistence (dark grey, bordered by red line) defines the region where the niche difference is greater than the fitness differences between species, implying that either species could invade from rare. The pale grey area denotes the region of 'priority effects,' where neither species could invade the other. b) Full posterior distribution of niche and fitness differences of one example species pair (D. sulfurigaster and D. simulans) using the best-fitting model. Note the irregular, non-circular, shape ultimately derived from correlations between fitted vales in the posterior distribution and alpha terms contributing to both niche and fitness differences. c) Equivalent contour plots for all 15 species combinations. Table 1 details the support assigned to each region for each pair under each treatment. Note that because the best model did not include varying competitive coefficients (α) between treatments, only the fitness difference (y-axis) shifts between treatments. An equivalent diagram using a model that allows α to vary is shown in Fig. S4.15.  Table 1. Summary of the posterior distribution across the coexistence plane between the two treatments for each species pair. The percentage of posterior draws that fall into the four major regions of the coexistence can provide a summary of the confidence in the outcome of a particular competition, given the data. The distinction between species 1 and 2 is arbitrary.

Discussion
Our results demonstrate the power of using a Bayesian approach within the modern coexistence framework, and in this case an inconsistent impact of consumers on coexistence. Natural enemies in the form of parasitoids affected species growth rates, but had no identifiable impact on pairwise competitive coefficients. Across our six species, no general susceptibility -fecundity trade-off was clearly observable and the presence of parasitoids did not lead to a consistent shift towards (or away from) a state of potential coexistence. In some cases, counter-intuitive results were observed. For instance, D. simulans and D. pandora show a negative relationship between fecundity and susceptibility (Fig. 2), but parasitism increases the likelihood of coexistence (Fig. 3c). This work highlights the importance and challenges of accurately measuring competition terms (α) to understand the consequences of susceptibility -fecundity trade-offs for coexistence (Petry et al. 2018).

Advantages of a Bayesian Approach
Fitting parametric models to ecological data is almost always both challenging and illuminating in equal measure. The framework of modern coexistence theory provides particular difficulties, since the quantities of interest are highly entangled (Barabás et al. 2018;Song et al. 2019). Because the fitted values of the parameters inevitably correlate to some degree, propagating parameter uncertainty from a multi-parameter posterior distribution gives a better understanding of the uncertainty associated with the key quantities. Such outputs are inherently multi-dimensional, and can pose considerable challenges both to intuition and to visualise. By contrast, previous approaches that propagate standard errors associated with individual parameters (e.g. Matías et al. 2018) or estimate uncertainty by bootstrapping observations (e.g. García-Callejas et al. 2020) are unable to capture these effects, and risk misrepresenting the confidence that can be placed in results.
Although we took an explicitly Bayesian approach, similar methods could be followed using multidimensional likelihood surfaces computed with a frequentist approach.
Understanding how the interrelationships between parameters affect uncertainty in the derived quantities based upon them is difficult a priori, since the contribution of different correlations will be contingent on both the underlying ecology and the experimental design. In the context of modern coexistence theory, even a moderate confounding of underlying growth rates and intra-specific competition terms may have a significant effect on the conditions for coexistence. In our study, the boundary line demarking the region of coexistence tended to be somewhat 'repelling'. This increases the confidence with which pairs of species are located within regions of the coexistence plane. Although our individual estimates of model parameters may be rather uncertain, we can often be comparatively confident in the result of competition.
This repelling effect can be attributed to indeterminacy between growth rates and intra-specific competition terms. One consequence of this can be seen clearly when the condition for coexistence (equation 6) is re-arranged in a simplified form as: Where there are strong positive correlations in the posterior draws of α and (shown in bold), this would tend to affect both sides of the inequality simultaneously, increasing the tendency that it maintains its value. However, it is important to note that there are a number of other correlations between core model parameters. For this reason, rather than trying to second-guess how these relationships may play out in other studies, we would recommend as far as possible propagating all uncertainties through to the final reported result.
Modern coexistence theory defines a sharp line between coexistence and the inability to invade from rare. This mathematical clarity and precision represents a major advantage, but poses challenges to the incorporation of inevitable experimental uncertainty. It is important to reiterate that the uncertainty of coexistence we discuss should be interpreted as statements of uncertainty in the underlying data, rather than implying a probabilistic outcome. The precision of estimates of parameters can differ between treatments, which can have the effect of reducing the support for coexistence, even if the median estimate does not change. Furthermore, these values should not be directly interpreted as some measure of 'distance from coexistence', although the concepts may be related. The framework of structural stability is better equipped to address such questions (Saavedra et al. 2017;Godoy et al. 2018).

Impact of natural enemies on competition coefficients
In natural systems, both direct effects of mortality and impacts on the interaction between species will determine the impact of consumers on species coexistence (Pringle et al. 2019). Our study included species that span a range of body sizes (0.94 mg -2.64 mg, mean wet mass of female adults) but also included two pairs that are very closely related phylogenetically (SI 1, Table S1.1). We expected that these choices would structure the competition outcomes, but we did not observe any clear patterns in the competition coefficients. Our study design allowed parasitoids to travel between vials, opening the prospect of either associational susceptibility or resistance mediated through parasitoid foraging behaviour (Holt & Kotler 1987;Abrams 2010). If competitive interaction modification effects were strong enough to be detectable, these would manifest as shifts between treatments in the fitted α values, and result in 'side-to-side' movement in the coexistence plane (e.g. SI 4, Fig. S4.15). However, our model selection procedure suggested that modelling changes to the competition coefficients matrix between treatments leads to overfitting. Nevertheless, our model selection test is a rather blunt all-or-nothing approach. It may well be that only particular competitive interactions are strongly affected by the parasitoid. Indeed, in our full model, while most of the between-treatment α differences cluster around zero, several show large divergences (SI 4 Fig. S4.13). There is no clear predictor for the terms that show these large shifts, occurring between different species. A more complex model fitting approach incorporating regularisation could potentially determine whether such cases are identifiable effects rather than model fitting noise, but is beyond the scope of this paper. Mechanisms that can be modelled as shifts in single species process (e.g. ) are inherently more identifiable statistically than changes to multi-species processes (e.g. α). In our models, this can be seen in that uncertainty in fitness differences was driven largely by the competitive response ratio √ , rather than the demographic ratio −1 −1 (SI 4, Fig. S4.3).
Density-mediated apparent competition derived from an increased total population of natural enemies, is not directly detectable with our single-generation design. If we had run a multigenerational experiment, the results could be expected to be complex: Sommers & Chesson (2019) show that the relative balance between apparent and resource competition is influenced by prey defensive actions. Methods to directly compare resource competition and predation effects in parameterised models are developing (Shoemaker et al. 2020) and will be an exciting area for future research.

Limitations of simple laboratory experiments
The modern coexistence theory approach is agnostic with respect to the particular factors involved, but will always be limited by the scope of the empirical data available. Experiments such as this one can rarely replicate full natural lifecycles of animals or capture the full range of conditions that occur throughout their habitat (Hart et al. 2018); and impacts on one life stage can mitigate differences at another (Moll & Brown 2008). Our methods likely favour species with a high egg-laying rate over those with a long reproductive lifespan. Equally, our experiment uses a highly simplified representation of the consumer pressure experienced by natural communities of Drosophila. We included a single pupal parasitoid species; in their natural habitats these species exist in a complex trophic network (Jeffs et al. 2020) including specialist parasitoids and larval parasitoids (which may directly influence larval competitive performance). Defences against parasitism may be specific to particular enemies, or behaviour-based and not observable in the simplified environment of the fly vials.
The simplicity of the experimental set-up cannot capture other plausible mechanisms that could be hypothesised to aid coexistence, such as environmental stress resistance (Chesson 2000;Kuang & Chesson 2010). It necessarily excludes a number of likely-significant environmental heterogeneities. For example, D. simulans is a cosmopolitan species that is not native to the Australian rainforest. It is strongly associated with human habitation (Parsons & Bock 1979), and may struggle to compete in natural forests against competitor species. Coexistence may only be possible at particular scales (Hart et al. 2017), or through spatial mechanisms not reflected in our experiment. In particular, we do not include any sense of dispersal capacity, or reproductive lifespan. Ultimately, our experiment can only represent a single snapshot combination of environmental conditions. Despite these limitations, we believe that our experiment will still be instructive in isolating particular aspects for more further investigation.

Scaling up from pairwise to community coexistence
The Drosophila species we investigated co-occur in their natural habitat of Australian rainforest (Jeffs et al. 2020), although it is important to note that this is not necessarily evidence that they can coexist in the technical sense (Siepielski & McPeek 2010). Furthermore, pairwise coexistence is not the same as community-level coexistence, and coexistence or otherwise between pairs does not necessarily translate to the multispecies community (Song et al. 2019), although they are likely to be at least somewhat related (Friedman et al. 2017;Broekman et al. 2019). A diversity of mechanisms can allow multispecies coexistence where pairwise coexistence is not possible Letten & Stouffer 2019; McPeek 2019).
The general approach of propagating joint-posteriors into derived metrics of coexistence from model terms could equally apply to other coexistence modelling frameworks (Godwin et al. 2020;Spaak & de Laender 2020) or alternative mechanisms that can facilitate community-level coexistence, such as intransitive competition (non-hierarchical competition, akin to an extended rock-paper-scissors tournament). We demonstrate this in SI 6 with a simple analysis of the impact of natural enemies on intransitivity. We did not identify a large amount of intransitivity, and it was not changed consistently by the presence of parasitoids. Although not a main focus of our investigation, our results add to a growing body of work suggesting limited impacts of intransitivity on community coexistence ).

Conclusion
Our experiment demonstrates that natural enemies can change the result of competition between species pairs, but that this effect is not consistent or necessarily predictable from pairwise fecunditysusceptibility relationships, reiterating the value of analysing parametrised models (Siepielski & McPeek 2010). Further work will be needed to investigate whether trait differences impact the niche and fitness differences observed in our system, and more experimental studies of the impact of natural enemies on coexistence are needed to form the basis for a more general understanding of their impacts. Modern coexistence theory is emerging as a powerful tool to investigate a range of biotic and abiotic impacts on coexistence experimentally (Lanuza et al. 2018;Matías et al. 2018), However, since the precision with which growth rates and competition coefficients can be estimated is limited (even within large studies such as ours), the fundamental inseparability of niche and fitness differences can impact the likelihood assigned to coexistence between any given species pair.
Although these challenges are significant, we hope our approach demonstrates that they can be readily and transparently addressed.