Natural enemies have inconsistent impacts on the coexistence of competing species

1. 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 the impact of natural enemies. 2. We tested experimentally the capacity for a generalist enemy to promote coexistence of competing insect species, and the extent to which any impact can be predicted by trade- offs between reproductive rate and susceptibility to natural enemies. 3. We used experimental mesocosms to conduct a fully factorial pairwise competition experiment for six rainforest Drosophila species, with and without a generalist pupal parasitoid. We then parameterised models of competition and examined the coexistence of each pair of Drosophila species within the framework of modern coexistence

2. We tested experimentally the capacity for a generalist enemy to promote coexistence of competing insect species, and the extent to which any impact can be predicted by trade-offs between reproductive rate and susceptibility to natural enemies.
3. We used experimental mesocosms to conduct a fully factorial pairwise competition experiment for six rainforest Drosophila species, with and without a generalist pupal parasitoid. We then parameterised models of competition and examined the coexistence of each pair of Drosophila species within the framework of modern coexistence theory. 4. We found idiosyncratic impacts of parasitism on pairwise coexistence, mediated through changes in fitness differences, not niche differences. There was no evidence of an overall reproductive rate-susceptibility trade-off. Pairwise reproductive rate-susceptibility relationships were not useful shortcuts for predicting the impact of parasitism on coexistence. 5. Our results exemplify the value of modern coexistence theory in multi-trophic contexts and the importance of contextualising the impact of generalist natural enemies to determine their impact. In the set of species investigated, competition was affected by the higher trophic level, but the overall impact on coexistence cannot be easily predicted just from knowledge of relative susceptibility.
Methodologically, our 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. This method should be widely applicable for understanding species coexistence in a range of systems.

| INTRODUC TI ON
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 species that would otherwise be competitively dominant (Paine, 1966), but their impact can be highly variable (Chase et al., 2002).
Analyses investigating 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; however, 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 it 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 modern coexistence theory, interactions mediated through competition for resources and those mediated by natural enemies in the form of 'apparent competition' are conceptually equivalent (Chesson, 2018;Chesson & Kuang, 2008), 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 interaction terms. Alongside the theory, empirical studies using the modern coexistence framework to investigate the impact of consumers on coexistence are now emerging .
Natural enemies can affect both fitness differences and niche differences between species. A frequently invoked mechanism is that they reduce fitness differences, mediated by a trade-off between population growth rate and susceptibility to attack (Kraaijeveld & Godfray, 1997;Paine, 1966). 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). Natural enemies can also potentially influence niche differences through modifying the overall competitive pressure exerted between competitors. Enemies may cause competitors to change their behaviour, or natural enemies may change their foraging behaviour as the relative abundances of the competing species change (Abrams, 2010).
Experimental studies of insect communities have considerable potential to further our understanding of the role of natural enemies in species coexistence. 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). However, to date, most experimental work within a formal MCT framework has been conducted on plant or microbial systems , including investigation of the effects of seed predation (Nottebrock et al., 2017;Petry et al., 2018) and seed pathogens (Mordecai, 2013). Although ecologists have been competing insect species (especially Drosophila) against each other for decades (Davis et al., 1998;Gilpin et al., 1986;Pearl, 1932;Worthen, 1989), very few studies have explicitly tested for mutual invasibility in insect systems (Siepielski et al., 2018), rather than particular conditions necessary for coexistence (Godwin et al., 2020;Spaak & de Laender, 2020).
Here we present an analysis parameterising the impact of generalist 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. We examine how niche and fitness differences are influenced by the presence of the parasitoid with a Bayesian approach that propagates uncertainty in model parameters through to niche and fitness differences. We show that the generalist natural enemy can influence the fitness differences between species but has inconsistent impacts on pairwise coexistence that would not be identifiable solely from reproductive ratesusceptibility relationships.

| Study species
We investigated a system of six naturally co-occurring Drosophila (SUL)) from a well-characterised Australian rainforest community (Hangartner et al., 2015;Jeffs et al., 2021;O'Brien et al., 2017), along with a co-occurring generalist parasitoid wasp from the genus Trichopria (Hymenoptera: Diapriidae, lab strain 66LD, for taxonomic details see Supporting Information 1 and Lue et al. (2021)). This parasitoid attacks Drosophila species at their immobilised pupal stage.
Successful development of a parasitoid results in the death of its host. Drosophila can mount a variety of physiological defences to parasitoid attack, which can be costly (Fellowes & Godfray, 2000) and similar defence-reproduction trade-offs have been identified in other insects (Schwenke et al., 2016).
Drosophila species have long been used in experimental ecology since they are straightforward to culture under laboratory K E Y W O R D S Bayesian, competition, Drosophila, model fitting, modern coexistence theory, natural enemies, parasitoid 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 northern Queensland, Australia (Jeffs et al. 2021). The six species span a range of body sizes (0.94-2.64 mg, mean wet mass of female adults) and phylogenetic divergences (Supporting Information 1, 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., 2021).

| Experimental design
We conducted a large, single-generation pairwise competition assay using a response-surface design following the guidelines of Hart et al. (2018), and in the presence or absence of pupal parasitoids. We

| 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. On the day preceding the start of the experiment, the founding flies were lightly anesthetised under CO 2 to be separated by sex. 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 under light CO 2 anaesthesia 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 hr before being removed. Figure 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 (35 cm × 20 cm × 14 cm) containing up to 88 vials (Supporting Information 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. This ratio of parasitoids to Drosophila was selected based on prior experience to result in a parasitism rate in line with field-parasitism levels of around 15% or greater (Jeffs et al. 2021).
The wasps were free to move between vials within each box. While it is possible that late-pupating Drosophila larvae could cross from one F I G U R E 1 Summary of the experimental regime showing the founding process, treatment and data collection procedures described in the main text 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 developed on the outside of the vials, suggesting that movement was negligible. After 30 hr, 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 3 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 (Figure 1). There was a highly consist-

| 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). We take as our underlying model of growth rates the Beverton-Holt model, without any carry-over between generations: where N i,t and N i,t+1 are the size of the population of species i in the founding and emerging generations, respectively. R 0i is a speciesspecific generational reproduction term, ii represents the intraspecific competitive effect of species i and ij are competition coefficients terms representing the effect of species j on i. It is important to emphasise that these α terms capture all impacts between species, including those through indirect means. The same function has been used for annual plants that have been the focus of most recent coexistence work Levine & HilleRisLambers, 2009;Pérez-Ramos et al., 2019).
The mean number of adults emerging per founder was consistently lower in vials with a solitary founding pair compared to those with three founders (Supporting Information 2, Figure S2.2).
Drosophila use cues from other females when determining oviposition sites (Duménil et al., 2016;Wertheim et al., 2002), which could be causing this Allee effect. It is not possible to incorporate such positive density dependence effects directly 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 excluded data from all singleton pair trials from the model fitting.

| 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). All code and data are archived . We assumed a negative binomial error distribution and fitted separate overdispersion terms for each parasitism treatment. Since we counted both species, each two-species experiment contributed two data points. Although this contributes a degree of non-independence, the gain in total sample size is significant and there was no discernible within-vial correlation in residual error (Supporting Information 2, Figure S2.3). Since we excluded the single-founder vials, our effective total 'n' across both treatments was 1,476, to estimate a maximum of 86 free parameters in our most complex model (including overdispersion terms).
The competition α terms were each fit with an underlying Gaussian prior (mean = 0, σ = 1) and constrained to be positive, implying only competitive interspecific interactions, with no facilitation possible. Growth rate terms R 0 were constrained to be positive, and fit with a weak Gaussian prior (mean = 10, σ = 10). Tests fitting α values with different constraints and priors did not give notably different results (Supporting Information 3).
We fitted and compared a sequence of candidate models of increasing simplicity. In the first, separate R 0 and α terms were fitted for each treatment. In the second, a single set of α terms, but separate R 0 terms were fitted across both treatments (i.e. parasitism affects emergence rate, but not competition coefficients). In the third model, a single set of parameters were fitted 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 Supporting Information 4.

| Pairwise coexistence criteria
The invasibility criterion for species coexistence requires each species to have a positive long-term low-density growth rate in the presence of the rest of the community (Chesson, 2000). Following Godoy and Levine (2014), coexistence of a pair of species is determined to be possible, following Beverton-Holt model dynamics,

when:
Equation 2 expresses the coexistence criterion as a balance between niche overlap (ρ) and fitness differences ( i j ) where: and Following Chesson (2000) and Godoy and Levine (2014), we define the niche difference between two species as 1 − . Certain α values can result in a niche overlap > 1, implying a negative potential for stabilisation through niche differences and highlighting possible cases of priority effects (Ke & Letten, 2018). Although some previous authors have elected to define such niche differences as 0, we present the values as calculated. Under either approach, coexistence is not expected under the MCT framework.
The Drosophila species we examined vary in their time to maturation. However, it is important to note that the mathematical expression of the invasibility criteria in a constant environment, Equation 2, is not dependent on generation time. That is to say, the sole criterion for the invader species to expand its population while the resident is at an equilibrium is that its population size in its next generation will be larger, regardless of how long a generation takes.
For each pair of species, we calculated the niche and fitness differences across 3,000 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, exclusion of either species or a situation of priority effects (see Figure 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.

| RE SULTS
Based on model comparison, there was strong support for the expectation that parasitism reduces the growth rate of each species, but little support for differences in competition coefficients between parasitism and no-parasitism treatments. Our second model   Figure 2b), the median value did not cross the boundary.
Strong positive correlations between fitted parameters, led by relationships between R 0i and ii terms, resulted in non-independence of estimates of the fitness and niche differences and non-elliptical posterior distributions (Figure 2b,c).
The parasitism treatment reduced population growth rates, but we did not identify a reproduction rate-susceptibility tradeoff across our six species (Figure 3). The pairwise relationship between reproductive rate and susceptibility to parasitism (proportional reduction in R 0 ) did not provide a useful heuristic in identifying the impact of parasitism on reducing fitness differences. Positive and negative reproduction rate-susceptibility relationships were both associated with increases and decreases in fitness differences.

| D ISCUSS I ON
Our results demonstrate an inconsistent impact of natural enemies on the pairwise coexistence of a guild of competing species.
Drosophila species showed differential susceptibility to parasitism, which manifests in shifts of the fitness differences in the coexistence plane. Natural enemies shift the fitness difference in favour of the species that is less susceptible, but this does not necessarily enhance coexistence. Across our six species, no general reproductive rate-susceptibility trade-off was observable and the presence of (2) Fully quantifying pairwise competitive ability requires considerable effort (Hart et al., 2018) 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 competition coefficients (α) between treatments, only the fitness difference (y-axis) shifts between treatments TA B L E 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 plane can provide a summary of the confidence in the outcome of a particular competition, given the data coexistence when the more susceptible species is otherwise competitively dominant, not just more fecund. This underscores the importance and challenges of accurately measuring competition terms (α) to understand the consequences of any given factor's influence on coexistence.

| Impact of 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 design allowed parasitoids to travel between vials, opening the prospect of either associational susceptibility or resistance mediated through parasitoid foraging behaviour (Abrams, 2010;Holt & Kotler, 1987). If these competitive interaction modification effects were strong enough to be detectable, they would manifest as shifts in the fitted α values and change the niche differences between parasitism and no-parasitism treatments. However, our model selection procedure suggested that modelling changes to the competition coefficients matrix between treatments leads to overfitting, that is, there was no identifiable impact of parasitism on pairwise competition coefficients. Resolving the relative importance of direct effects of predators versus their indirect effects on other species is a long-standing challenge (Werner & Peacor, 2003). Our study suggests that in this case the direct effects are dominant. However, mechanisms that can be modelled as shifts in single-species parameters (e.g. R 0 ) are inherently more identifiable statistically than changes to multi-species processes (e.g. ij ) as they can be determined more directly from simpler experiments.
Nevertheless, our model selection test is a rather blunt all-ornothing approach. It may well be that only particular competition coefficients are strongly affected by the parasitoid. Indeed, in our Model 1, while most of the between-treatment α differences cluster tightly together, several show large divergences (Supporting Information 4, Figure S4.14). 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.
Our single-generation design meant that density-mediated apparent competition derived from an increased total population of natural enemies was not directly detectable. If we had run a multigenerational experiment, the results could be expected to be complex: Sommers and Chesson (2019) show that the relative balance between apparent and resource competition is influenced by prey defensive actions. Methods to compare resource competition and predation effects directly in parameterised models are being developed (Shoemaker et al., 2020) and will be an exciting area for future research.

| Model fitting considerations within coexistence theory
Recent meta-analyses have emphasised the importance of parameterised models for understanding mechanisms of coexistence (Broekman et al., 2019), particularly since large-scale global environmental changes will create novel communities (Alexander et al., 2016). 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, even before comparison to field contexts. 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.  Saavedra et al., 2017).
Nonetheless, 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, but 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.
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 intraspecific competition terms may have a significant effect on the conditions for coexistence. In our study, the posteriors were very often curved-posterior draws generating large niche differences were associated with larger relative fitness differences. This has the beneficial effect of often increasing 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 effect can be attributed to indeterminacy, and hence correlation in the posterior, between growth rates and intraspecific competition terms. One consequence of this can be seen clearly when the condition for coexistence (Equation 2) is re-arranged in a simplified form as: In our data, there are strong positive correlations in the posterior draws of ii and R 0i (highlighted). Hence, uncertainty would tend to affect both sides of the inequality simultaneously, increasing the tendency that it maintains its value.

| 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 life cycles of animals or capture the full range of conditions that occur throughout their habitat (Hart et al., 2018). Impacts on one life stage can mitigate differences at another (Moll & Brown, 2008) and our method does not capture differences in reproductive life span. 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., 2021) including specialist parasitoids and larval parasitoids (which may directly influence larval competitive performance).
Coexistence may only be possible at particular spatial scales (Hart et al., 2017), or through spatial mechanisms not reflected in our experiment. The laboratory experiments cannot capture environmental heterogeneities hypothesised to aid coexistence (Chesson, 2000;Kuang & Chesson, 2010). Despite these limitations, we believe that our experiment is still instructive in isolating particular aspects for further investigation.
The Drosophila species we investigated co-occur in their natural habitat of Australian rainforest (Jeffs et al., 2021), although it is important to note that this is not necessarily evidence that they can coexist in the sense of mutual invasibility (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 multi-species community (Song et al., 2019), although they are likely to be at least somewhat related (Broekman et al., 2019;Friedman et al., 2017). A diversity of mechanisms can allow multi-species coexistence where pairwise coexistence is not possible (Letten & Stouffer, 2019;Levine et al., 2017;McPeek, 2019).

| CON CLUS IONS
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 reproduction rate-susceptibility relationships, reiterating the value of (5) 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. 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 beyond plants (Viola et al., 2010), but we did not find evidence to support the concept that generalist enemies mediate community coexistence. 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.

ACK N OWLED G EM ENTS
We thank the Hrček lab and the Oxford Fly group for advice and use of facilities, Chia-Hua Lue for advice on parasitoid taxonomy, František Sládeček and Jamal Holle for assistance maintaining cultures, and Mark Wong for constructive comments on the manuscript.
We also thank the anonymous reviewers for highly constructive comments. Drosophila and parasitoid cultures were exported under permit PWS2016-AU-002018 granted by the Department of Energy and Environment of the Australian Government.