Connecting higher‐order interactions with ecological stability in experimental aquatic food webs

Abstract Community ecology is built on theories that represent the strength of interactions between species as pairwise links. Higher‐order interactions (HOIs) occur when a species changes the pairwise interaction between a focal pair. Recent theoretical work has highlighted the stabilizing role of HOIs for large, simulated communities, yet it remains unclear how important higher‐order effects are in real communities. Here, we used experimental communities of aquatic protists to examine the relationship between HOIs and stability (as measured by the persistence of a species in a community). We cultured a focal pair of consumers in the presence of additional competitors and a predator and collected time series data of their abundances. We then fitted competition models with and without HOIs to measure interaction strength between the focal pair across different community compositions. We used survival analysis to measure the persistence of individual species. We found evidence that additional species positively affected persistence of the focal species and that HOIs were present in most of our communities. However, persistence was only linked to HOIs for one of the focal species. Our results vindicate community ecology theory positing that species interactions may deviate from assumptions of pairwise interactions, opening avenues to consider possible consequences for coexistence and stability.

Lotka-Volterra competition or predator-prey models assume that the direct interaction between two species is an intrinsic property, thus the strength and sign of the interaction are independent of the community in which these species are embedded (Werner & Peacor, 2003).However, when more than two species co-occur (Worthen & Moore, 1991), indirect effects can arise through chains of pairwise interactions, or through higher-order interactions (HOIs), that is, changes of the per capita effect of one competitor on another in the presence of additional species (Levine et al., 2017;Wootton, 1994).Since interaction chains and HOIs emerge from fundamentally different mechanisms, it is important to understand their underlying mechanisms for appropriate modeling and inference (Levine et al., 2017).
Interaction chains emerge when pairwise interactions are embedded in a network of interactions.A series of such direct interactions between species pairs can lead to connections between species that do not directly interact with each other (Wootton, 1993).Despite introducing indirect pathways between species, interaction chains are the result of a series of fixed, strictly pairwise interactions.In contrast, "higher-order interactions" (HOIs (Case & Bender, 1981;Letten & Stouffer, 2019;Mayfield & Stouffer, 2017); also called "interaction modifications" (Abrams, 1983)), encompass the nonadditive effects of a species on the pairwise interaction between a focal pair.Conventionally, the species altering the interaction is a third species (i.e., an interspecific HOI); however, it can also be one of the focal species (i.e., an intraspecific HOI) sensu Letten and Stouffer (2019).Wootton (1994) observed both types of indirect effects in the upper zone of a rocky intertidal community.As an example of an interaction chain, bird predators indirectly increase acorn barnacle (Balanus glandula) abundance by consuming limpets (Lottia digitalis), which dislodge or consume young acorn barnacles.An example of an interaction modification, or interspecific HOI, is barnacles altering the bird-limpet interaction by changing the ability of birds to find limpets due to the similar color of L. digitalis and barnacles shells.HOIs are most commonly described in terms of how a third species impacts the interaction between a focal species pair.However, intraspecific HOIs are also important to consider and represent the cumulative impacts of interactions among the individuals of each "competitor species" (that is, intraspecific crowding) on the focal species (Letten & Stouffer, 2019;Mayfield & Stouffer, 2017).In the intertidal community, intraspecific HOIs could occur when birds influence the foraging behavior of barnacles, causing them to aggregate in areas sheltered from bird predation, which leads to increased intraspecific competition.Interaction chains can be predicted with only a knowledge of pairwise species interactions; in contrast, descriptions of intra-and interspecific HOIs require a knowledge of all species combinations involved (Wootton, 1994).
Theoretical work has suggested that empirical studies of interactions between species should consider HOIs (Grilli et al., 2017), given that simple models of pairwise interactions fail to explain the stable persistence in simulation models of very large ecological communities (Barabás et al., 2016;Clark, 2010;Gibbs et al., 2022;Kleinhesselink et al., 2022;Levine et al., 2017).Recently, many theoretical advances have been made to understand the influence of HOIs on community dynamics (Gibbs et al., 2022;Grilli et al., 2017;Letten & Stouffer, 2019;Mayfield & Stouffer, 2017;Singh & Baruah, 2021).A particular focus has been on the effect of HOIs on the stability of communities, specifically the persistence of species in a community.However, simulation-based studies often consider highly complex and diverse communities (Bairey et al., 2016;Singh & Baruah, 2021), making it impractical to verify their findings through experimental manipulations.
A valuable approach to simplify the complexity of communities is to analyze their component community modules (Holt, 1997).
Studies of community modules have revealed how interactions between species affect species persistence (Kondoh, 2008;Mayfield & Stouffer, 2017;McCann et al., 1998).Theoretical investigations on food web modules have highlighted the crucial role of weak interactions, especially omnivorous links, in maintaining community stability (Emmerson & Yearsley, 2004;McCann et al., 1998).Experimental manipulations of these modules have confirmed that weak interactions, coupled with strong interactions mediated by generalist consumers, enhanced community stability by reducing interaction strength (Rip et al., 2010).Experimental work has indicated that HOIs weaken intense pairwise competition and predator-prey interactions, thereby stabilizing the community and promoting coexistence (Kratina et al., 2007;Sundarraman et al., 2020).
Despite the theoretical support for the importance of HOIs for coexistence and stability, empirical evidence connecting HOIs with ecological stability remains very scarce (but see Mayfield & Stouffer, 2017;Mickalide & Kuehn, 2019;Sundarraman et al., 2020).This is partly due to the challenges of empirically quantifying HOIs (Billick & Case, 1994), which requires the evaluation of the interaction strength between two species accounting for changes in density or the presence of additional species (Adler et al., 2018).Even in simple communities, this can quickly become logistically infeasible.
In addition, empirical research also requires detailed time series data to evaluate species persistence over extended periods of time, presenting a second logistical hurdle.

Our goal was to experimentally test the causal links between
HOIs and the persistence of species within communities.Microcosms are a convenient tool to investigate concepts in community ecology because they are easily manipulated, can be highly replicated, and have large population sizes (Altermatt et al., 2015).
Moreover, a range of analytical approaches is available to estimate intra-and interspecific interaction strengths within these systems (Carrara et al., 2015).We, therefore, used microcosms to experimentally examine HOIs and community persistence in aquatic microbial communities.We used six different community compositions, which all included a focal pair of consumers competing for the same resource.We tested the effect of two additional competitors as well as the presence of a predator on the interaction strength between the focal pair.Time series of population dynamics were collected to quantify the strength of species interactions and examine the effect of HOIs on persistence and address the following questions: 1. Are HOIs detectable in our communities?2. Do HOIs differ depending on the identity and trophic role of the third species?3. Is there a relationship between the presence of HOIs and population persistence?
We hypothesized that additional competitors would influence the effect of the two focal species on each other since all species engage in resource competition.For instance, in the presence of a third species, focal species may shift their foraging behavior and thereby their pairwise effect on one another.We also expected that higher-order effects would arise in communities with a predator whose behavior is affected by different consumer body sizes which could render one prey susceptible to predation while the other may be attacked but is too large to be consumed.All else being equal, we expected additional competitors to strengthen the competition for shared resources shortening population persistence, but additional predators to weaken competition between the focal pair of species, potentially extending species persistence.

| Food web construction and culture conditions
To experimentally test HOIs, we constructed an experimental food web.The ciliates of our microbial food web were primarily bacterivorous (Pennekamp et al., 2018); thus, the experimental communities were sustained on a mix of three species of bacteria (Bacillus subtilis, Serratia fonticola, and Brevibacillus brevis) decomposing the protozoan pellet medium (PPM; provided by Carolina Biological Supplies; concentration of 0.55 g/L, see Altermatt et al., 2015).(Woodruff & Spencer, 1921).
Due to the large body size of P. caudatum (170-300 μm) (Foissner & Berger, 1996) and S. teres (150-400 μm) (Bick & World Health Organization, 1972), we expected these species to interfere with the predator because they are more difficult or impossible to consume; however, depending on the mode of attack of the predator, even failed consumption could result in death of prey.
We chose these five eukaryotic ciliates (hereafter referred to by the genus name only) because their pairwise trophic interactions are well documented by previous studies (Daugaard et al., 2019;Pennekamp et al., 2018;Tabi et al., 2019).The consumptive and competitive interactions present in our experimental food web are shown in Figure 1.

| Experimental design
To assess how the interaction strength between two focal species is altered by adding a third or more species, we established the following treatments: (1) two prey (Colpidium and Dexiostoma) alone to investigate the direct effects between the two focal species (CD) (2) two prey and Paramecium (CDP).(3) two prey and Spirostomum (CDS).These communities were also grown in the presence of the top predator Spathidium (CDPd, CDPPd, CDSPd; predator = Pd) (Figure 2), yielding six treatments in total.For each treatment, we cultured four replicates resulting in 24 microcosms in total.
To start the experiment, we took out 20 or 30 mL of the bacterized PPM medium (depending on experimental community composition) and replaced it with 10 mL of the stock culture of each consumer species at their carrying capacity.Therefore, the starting  densities of all consumers were 10% of their carrying capacities.
In the predation treatments, 4 days after the introduction of the consumer species, ten predator individuals were pipetted from the maintenance plates to the microcosms.To assure establishment in each microcosm, another ten individuals of predators were added after 8 days.The experiments were conducted at 15°C with no light, and suitable growing conditions for ciliates (Altermatt et al., 2015).

| Video sampling and species classification
Experimental units were sampled every second day for 53 days to capture time series of the dynamical changes in the abundance of ciliates.For each sampling event, the microcosm was gently agitated, a subsample of 250 μL was mounted onto a glass slide and covered with a glass lid.Three five-second videos (at 25 frames per second) were taken using 25× magnification on a stereomicroscope (Leica M205 C) mounted with a digital CMOS camera (Hamamatsu Orca Flash 4.0 C11440, Hamamatsu Photonics, Japan) with dark field illumination.We took three subsamples for each microcosm to get a precise estimate of abundance and processed the videos with the R package BEMOVI (version 1.0.2) (Pennekamp et al., 2015).
We took the mean of the three videos as our measure of density.If no individuals were detected, we assigned a zero.The volume lost from the microcosm by sampling during the course of the experiment was replaced with a fresh bacterized PPM medium.
For species classification, we trained a support vector machine (SVM) classifier on 350 to 400 randomly chosen and manually labeled trajectories of each species across the community compositions and time with the R package e1071 (Meyer et al., 2022).We also included a "noise" class in our classifier representing spurious trajectories due to background movement.Twenty morphological and movement features extracted from an established classification pipeline (Pennekamp et al., 2017) were selected to train the SVM to distinguish among classes based on information about body size and movement patterns (Table S1).As ciliate phenotypes may change over time (Pennekamp et al., 2017), we included "week number" to enhance the accuracy of the classifier.Further details about the classification can be found in the Appendix S1.

| Abundance
To understand the effect of community composition on the abundance of the species, we calculated the mean abundance over the entire 53-day experiment for each species in each of the six communities.Mean abundance of each of the five species was analyzed with two-way ANOVA (predictors: identity of the third species [levels: Paramecium or Spirostomum] and the presence of the predator [levels: present or absent]).

| Estimating species interaction strengths
We estimated the effect of species j on the population growth rate of species i within a regression framework, where the population growth rate of species i is regressed against the densities of species i and j to determine the intra-and interspecific effects on the population growth rate of species i (Pfister, 1995).The per capita population growth rate of each species i was calculated as (t + 1) − t where N i,t+1 is the abundance at time t + 1 and N i,t is abundance at time t.We used the gauseR package to calculate the per capita population growth rate (Mühlbauer et al., 2020).To explore how interaction strength and HOIs may differ between communities, we compared the fit of six models to the per capita population growth rate of the focal species (either Dexiostoma or Colpidium).First, we fitted the Lotka-Volterra model with no HOIs: For the simple additive model, λ i is the intrinsic population growth rate and α ij represents both the intra-and interspecific interaction coefficients.To test for intra-and interspecific HOIs, we then added interaction terms to the additive model LV, resulting in the interactive LV model (Letten & Stouffer, 2019).For the interactive LV model only including intraspecific HOIs, we included the interaction term β ijj for the effect of species j on i, where β ijj captures the cumulative impacts of intraspecific interactions on the focal species: The identity of species j may or may not be the same as species i.
Next, we fitted an interactive LV model that only includes interspecific HOIs, where β ijk captures the cumulative impacts of interspecific interactions on the focal species: Here, the identity of species k strictly excludes the focal species i.
Finally, we tested the fit of the interactive LV model including both intra-and interspecific HOIs (full HOI model): Population dynamics for each experimental community.Rows are six different treatments in our experiment.Each treatment is represented by the species it contains (C for Colpidium, D for Dexiostoma, P for Paramecium, S for Spirostomum and Pd for predator).The yaxes are plotted on a logarithmic scale.Results for each replicate are shown in columns.Most species in competition treatments (first three rows) survived to the last day of measurement but rapidly dropped below detection limits when housed with predators (last three rows).
To test whether the shape of the density dependence deviates from the linear form of the additive Lotka-Volterra model, we fitted the Ricker model with nonlinear density dependence.First, we fitted the additive model without HOIs: and interactive form only including interspecific HOIs (terms as previously defined): We used generalized linear models with the "identity" link function assuming Gaussian errors to estimate the coefficients of the additive and interactive LV model, considering all errors as measurement errors.For the Ricker model, we used a GLM with the "log" link function, assuming also Gaussian errors.For the fitted models, the intercept can be interpreted as the maximum intrinsic population growth rate, while the regression coefficients represent the per capita effect of species i on itself and the per capita effects of species j on species i. Interaction terms describe the mediating effects that the density of species j can have on the effect of species i (and vice versa) on the population growth rate of the focal species (Letten & Stouffer, 2019).We performed an AICc model selection adjusting for the small sample size when identifying the model that best explained the population growth rate of each focal species (Burnham & Anderson, 2002).AICc was preferred over the Bayesian information criterion since the models which we fit are approximate and we do not necessarily expect to include the true model in our set (Aho et al., 2014).We then calculated the difference (ΔAICc) between the model with the lowest AICc and all other models.Models with a ΔAICc larger than 2 are considered to be different; thus, the model with the fewest parameters and a ΔAICc less than 2 is considered the most parsimonious model.Since the presence of predators drove some prey to early extinction, we only estimated interaction strengths in communities without predators.

| Estimating species survival (persistence)
We measured the persistence of the focal species Colpidium and Dexiostoma, as well as the predator Spathidium as a proxy of the persistence of these species using nonparametric survival analysis.We estimated the time to event (i.e., extinction), accounting for right-censoring, that is, the possibility of missing future extinctions since we only recorded the abundance for 53 days.We tested the effect of community composition on the Kaplan-Meier estimate of the focal species with log-rank tests, followed by pairwise tests of community composition in case of a significant community wide effect with the survminer package (Kassambara et al., 2021).
All of the above analyses including the species classification were performed with the statistical computing environment R (R Core Team, 2022).

| Population dynamics over time
The population dynamics of the protists differed depending on the community context (Figure 2

| Higher-order interactions
Model selection revealed that HOIs influenced the per capita population growth rates of Dexiostoma and Colpidium in some communities (Table 1; Figures S1-S7).For Colpidium, the interactive LV model only including intraspecific HOIs was best supported in the

| Population persistence
In the absence of predators, Colpidium survival was contingent on the community composition (log-rank test p = .0066,n = 4; Figure 4a).
Without predators, adding Spirostomum increased Colpidium survival (pairwise log-rank test p = .02,n = 4; Figure 4a), whereas adding Paramecium did not.In contrast, Dexiostoma persisted in the absence of predators regardless of the community composition (log-rank test p = .37,n = 4; Figure 4c).In the presence of predators, Colpidium survival did not depend on the identity of the competitor species (Figure 4b).In contrast for Dexiostoma, the community composition did impact survival (log-rank test p = .024,n = 4), since adding Paramecium to the community increased Dexiostoma survival (Figure 4d).

| Patterns of population abundance
Overall, species' abundances and coexistence were influenced by community composition.Colpidium and Dexiostoma grew to lower densities when cultured with Paramecium, a sign of competition for shared resources.Colpidium abundance was not affected by the presence of Spirostomum.However, the mean abundance of Dexiostoma was lower in the presence of Spirostomum, suggesting competition between the two species, but to a lesser degree than with Paramecium.These patterns are in line with known phylogenetic relationships between these species, where Colpidium is most closely related to Dexiostoma, then to Paramecium and finally to Spirostomum (Violle et al., 2011).These relationships have been found to predict competitive exclusion due to the similarity in mouth size between species, which in turn defines the feeding niche of species (Violle et al., 2011).
Predation lowered the abundances of the two focal prey species, but this effect could be partially mitigated in the presence

TA B L E 1 Model selection table
showing the AICc values for the per capita population growth rate of the two focal species across the three community compositions.

F I G U R E 4
why some species may interfere with predator foraging.Dexiostoma similarly benefited from Paramecium when the predator was present, reaching on average higher abundance than when Colpidium and Spirostomum were present.
Surprisingly, both Paramecium and Spirostomum densities decreased in the presence of the predator, indicating that the predator had a negative effect on the two largest protist species.
Furthermore, the abundance of the predator was greater when Paramecium was present compared to when it was cultured with only the focal species pair, or the focal pair and Spirostomum.It thus appears that Paramecium was consumed in addition to the two focal prey and that while Spirostomum was attacked, it was potentially too large to be consumed.This observation could be explained by the mode of prey capture employed by the predator Spathidium.The cell mouth (cytostome) of Spathidium is furnished with a rod-like tip (toxicyst) to paralyze or kill other microorganisms for easier consumption (Fyda et al., 2005).But while the toxicyst make it possible to paralyze large prey, phagocytosis still requires that prey organisms can be engulfed.

| Evidence for higher-order interactions
Comparing additive and interactive LV models revealed the significant effect of intraspecific HOIs on the population growth rate for both Colpidium and Dexiostoma in the two species communities (Figure 5a).While our experiment was not designed to test the underlying mechanisms that give rise to intraspecific HOIs, one plausible explanation could be the presence of an "information cascade," wherein individual behaviors within the same species are regulated by the actions of others (Bikhchandani et al., 1992;Potts, 1984).Studies have shown that the behavior of neighboring fish influences the direction and speed of the school (Ioannou et al., 2011).If a fish senses a signal of danger and turns, it creates a pressure wave in the water, and other fish respond to this pressure wave by turning as well (Ioannou et al., 2011).Although ciliates do not have as elaborate sensory systems as fish, they are capable of sensing changes in local population density and adjusting their movement strategies accordingly (Fronhofer et al., 2015;Pennekamp et al., 2014).Since movement and foraging are often closely linked (Van Dyck & Baguette, 2005), it is possible that movement-related crowding effects could lead to the emergence of intraspecific HOIs.
Interspecific HOIs were only found on the population growth rate of Colpidium when Dexiostoma and Spirostomum were present (Figure 5b).The prevalence of intra-versus interspecific HOIs is noteworthy, since most studies have so far focused on interspecific HOIs.
Work on plant communities has shown that the prevalence of intraor interspecific HOIs can vary across focal species and that neither type is generally more or less important (Mayfield & Stouffer, 2017).
Since coexistence will be determined by the relative strengths of intra-and interspecific competition (Gibbs et al., 2022;Singh & Baruah, 2021), changes in intraspecific interaction strength alone can influence the persistence of communities.We therefore recommend that empirical studies investigating HOIs embrace a definition of HOIs that includes both intra-and interspecific effects (e.g., Mayfield & Stouffer, 2017).

| The nature of the detected HOIs
Some authors have argued that HOIs are not ecological processes in their own right but are instead emergent properties of phenomenological models (Letten & Stouffer, 2019).For example, a Lotka-Volterra competition model imposes a linear relationship between the density of one species and the growth of the other.
However, if this relationship is not linear in nature, the addition of HOI terms may improve the fit, but the model would not accurately describe the species' interactions.Properly accounting for nonlinear density dependence is an important prerequisite to estimating interaction strength (Hart et al., 2018) and avoiding erroneous conclusions about the presence of HOIs when nonlinearity is present (Kleinhesselink et al., 2022;Letten & Stouffer, 2019).
To investigate the potential for, and source of, nonlinear density dependence, we fit models that included either linear (i.e., additive LV) or nonlinear (i.e., additive Ricker) competition coefficients, as well as nonlinear intra-and interspecific HOI terms (i.e., the interactive models).We observed nonlinear density dependence in all communities, as either the Ricker or the interactive LV model  et al., 2007;Terry et al., 2017).

| Limitations of our work
While our study provides insights into the occurrence and implications of HOIs in aquatic microcosms, the small sample size used in our experiment limits the generalizability of our findings.Detecting interactions usually requires larger sample sizes, unless deviations from the additive model are very large (Burgess et al., 2022).Due to the small sample size of our experiment, we may have missed important but weaker interactive effects.In addition, we were only able to examine the interaction of a single focal pair across a few community compositions.To completely disentangle the drivers of community persistence would require to measure all interaction strengths between species, which was logistically impossible in our study.Since our study only investigated HOIs between a focal pair, it is possible that we missed additional interspecific HOIs between the focal pair and the predator.Future research with a larger sample size and a wider range of competitors and predators is warranted to further elucidate the complex nature of HOIs in ecological communities.

| CON CLUS IONS
The role of HOIs in community dynamics remains a key issue in community ecology (Gibbs et al., 2022;Levine et al., 2017).Currently, very few empirical studies measure HOIs in communities of real species and study how HOIs affect stability.Our study demonstrates that HOIs can play a role for species persistence and provides some support for the stabilizing effects of HOIs on competitive communities shown in theoretical studies (Grilli et al., 2017).Our study further pro-

ACK N OWLED G M ENTS
We thank two reviewers and the associate editor for constructive feedback on the manuscript.Chenyu Shen, Kimberley Lemmen and Frank Pennekamp received financial support from the Swiss National Science Foundation (grant 310030_197811 awarded to FP).We thank Uriah Daugaard for inspiring the ciliate illustrations in Figure 1.

CO N FLI C T O F I NTE R E S T S TATE M E NT
None.

O PE N R E S E A RCH BA D G E S
This article has earned Open Data and Open Materials badges.

Four
bacterivorous ciliates (Colpidium striatum, Dexiostoma campylum, Paramecium caudatum, and Spirostomum teres) constitute the intermediate consumers in our experimental microcosms.Both C. striatum (length: 50-100 μm) and D. campylum (35-90 μm) are consumed by the top predator Spathidium sp.(40-300 μm), which cannot survive only on bacteria as prey ).When Colpidium and Dexiostoma competed only with each other, we consistently observed the extinction of Colpidium toward the end of the experiment, while Dexiostoma persisted.The addition of Paramecium or Spirostomum led to a longer period of coexistence of Dexiostoma and Colpidium.The addition of the top predator Spathidium considerably shortened the persistence of all species.Thus, the addition of a third competitor increased the period of coexistence of the focal species pair, whereas the addition of the predator led to the rapid collapse of the whole community.The addition of the competitor Paramecium had a negative effect on the mean abundance of both Colpidium (b = −0.60,95% CI = −0.89 to −0.32, n = 4) and Dexiostoma (b = −0.81,95% CI = −1.0 to −0.59, n = 4, Figure 3).In contrast, adding the competitor Spirostomum only had a negative effect on the mean abundance of Dexistoma (b = −0.43,95% CI = −0.67 to −0.19, n = 4, Figure 3).The addition of the predator had a negative effect on the mean abundance of all competing protist species (Figure3).In communities composed of just the focal pair, the presence of the predator resulted in a similar reduction in mean abundance for both Dexiostoma:(b = −1.8,95% CI = −2.0 to −1.5, n = 4) and Colpidium (b = −1.9,95% CI = −2.1 to −1.6, n = 4, Figure3).In communities with a third competitor, the predator had a greater negative effect on the mean abundance of Spirostomum (b = −1.7,95% CI = −1.9 to −1.5, n = 4) than Paramecium (b = −1.5, 95% CI = −1.7 to −1.2, n = 4, Figure3).The addition of Spirostomum had no effect on the mean abundances of Dexiostoma, Colpidium and the predator compared to the focal pair cultured with just the predator.Intriguingly, the presence of both Paramecium and the predator had a positive effect on the density of both Colpidium (b = 0.83, 95% CI = 0.43 to 1.2, n = 4) and Dexiostoma (b = 0.96, 95% CI = 0.62 to 1.3, n = 4; Figure 3) in comparison to when the focal species were exposed to the predator in isolation.This translated into higher density for the predator Spathidium when Paramecium, Colpidium and Dexiostoma were present (b = 0.66, 95% CI = 0.48 to 0.84, n = 4; Figure 3).

F
Density of focal species as a function of community composition.Each panel shows one of the protist species.Whether the predator was present is shown on the x-axis.The treatments with predators are shown with triangles, and the treatments without predators are shown with circles.The three competitive communities are shown in color (CD = Colpidium & Dexiostoma, CDP = Colpidium, Dexiostoma & Paramecium, CDS = Colpidium, Dexiostoma & Spirostomum).The solid shapes show the mean, and the error bars show the standard error of the mean (calculated across the duration of the experiment), and the transparent shapes show the actual data points.
Spirostomum most likely changed the effect of Dexiostoma on Colpidium, because both Dexiostoma and Colpidium strongly competed for the same resources, whereas the direct effect of Spirostomum on Colpidium was more limited.It is possible that there is an overlap in the size of the bacteria consumed by Dexiostoma and Spirostomum, while Colpidium uses a different size class.In such a case, Dexiostoma would compete with both species but Colpidium would only compete with Dexiostoma.This would also explain why we did not observe an interspecific HOI on the population growth rate of Dexiostoma when Spirostomum and Colpidium were present.Interestingly, there was no HOI of Paramecium on either Dexiostoma or Colpidium, although Paramecium was a stronger competitor than Spirostomum.Previous work has found strong evidence that interspecific HOIs can affect the coexistence and dynamics of aquatic microbial communities (Mickalide & Kuehn, 2019).Mickalide and Kuehn (2019) observed that Escherichia coli can invade cultures of the alga Chlamydomonas reinhardtii or the ciliate Tetrahymena thermophila but fails to invade a community where both species are present.The invasion resistance of the algae-ciliate community arises from an HOI caused by the algal inhibition of bacterial aggregation, which leaves bacteria vulnerable to predation.

F
I G U R E 5 Higher-order interactions detected in our experiment.Direct pairwise interactions are shown as red arrows, HOIs are shown in blue.(a) In the two-species community of Colpidium and Dexiostoma, intraspecific HOIs were detected, indicating Dexiostoma modified the effect of Colpidium on itself while Colpidium modified the effect of Dexiostoma on itself.(b) In the community with Colpidium, Dexiostoma, and Spirostomum, an interspecific HOI was detected (the name of the affected species is marked with a light-colored box), suggesting Spirostomum modified the interaction between the two focal species.
had the lowest AICc values.The observation that models which included HOIs had the lowest AICc suggests that there are processes resulting in nonlinearity which the simple additive LV and the nonlinear Ricker model are missing.Therefore, we believe the HOIs detected are not an artifact of nonlinear density dependence, but true behavioral responses.A mechanistic model that allows for such nonlinearity may provide further insights into the ecological processes driving community dynamics.4.4 | Do HOIs affect the abundance and persistence in competitive communities?Persistence was different for the two focal prey species.Dexiostoma showed high persistence regardless of the additional species present, while Colpidium had higher persistence in the presence of a third species.The increased persistence of Colpidium when cultured with Paramecium was not driven by an HOI, since the additive Ricker model was best supported by the data in the CDP community.But the highest persistence of Colpidium was detected in the presence of Dexiostoma and Spirostomum, where the interactive model was best supported.This pattern is consistent with a stabilizing effect due to interspecific HOIs.Dexiostoma showed a similar persistence across all combinations of competitors, despite a nonsignificant trend to lower persistence when Paramecium was present.Dexiostoma also showed variation in the presence of HOIs, but the pattern is not consistent with a stabilizing role of HOIs: when present with only Colpidium, intraspecific HOIs were observed and Dexiostoma persisted till the end of the experiment.When Paramecium or Spirostomum was added, HOIs were no longer detected but neither did the persistence change.This suggests that changes in HOIs do not always result in changes in the persistence of focal species.The predator itself showed the highest persistence when feeding on communities that contained Colpidium, Dexiostoma and Paramecium.When feeding on the focal pair without competitors, it showed intermediate persistence times, while communities of the focal pair with Spirostomum persisted the least well.Hammill et al. (2015) found a similar effect for nontrophically interacting species in a food web: if species that do not interact trophically are present, prey persists for longer in diverse food webs.Since we could not quantify the interaction strength of the predator on Colpidium and Dexiostoma in the presence of Paramecium or Spirostomum, we could not determine whether the change in predator persistence was due to HOIs or simply the direct and indirect effects among the members of the food web.Quantifying the role of trophic interaction modifications with functional responses is an exciting area for future research(Kratina vides an example of how to identify HOIs by combining experiments with time series analyses and hence pave the way for more studies that study the occurrence and consequences of HOIs in natural communities.Our study adds to the growing body of research that shows that species interactions can deviate from the pairwise expectations when additional species are present, and a new theory needs to explore the consequences for community coexistence and stability.AUTH O R CO NTR I B UTI O N S Chenyu Shen: Conceptualization (equal); data curation (equal); formal analysis (equal); investigation (equal); methodology (equal); validation (supporting); visualization (equal); writingoriginal draft (equal); writing -review and editing (equal).Kimberley Lemmen: Conceptualization (supporting); formal analysis (supporting); investigation (equal); methodology (equal); writingoriginal draft (supporting); writing -review and editing (equal).Jake Alexander: Conceptualization (supporting); investigation (supporting); methodology (supporting); supervision (supporting); writing -original draft (supporting); writing -review and editing (supporting).Frank Pennekamp: Conceptualization (equal); data curation (equal); formal analysis (equal); funding acquisition (lead); investigation (equal); methodology (equal); project administration (lead); resources (lead); supervision (lead); validation (lead); visualization (equal); writing -original draft (equal); writing -review and editing (equal).
(Hansen et al., 1994)from the presence of Paramecium when predators were present, having a higher average abundance than would be expected if the effects of both predation and competition were additive.This might be explained by the interference of Paramecium with the predator, indirectly leading to a decrease in interaction strength between the prey and the predator.Ciliate predators show a preferred predator-prey size ratio of 8:1(Hansen et al., 1994), explaining Additive and interactive versions of the Lotka-Volterra (LV) and Ricker models were fitted to the data.K is the number of estimated parameters, ΔAICc is the difference between the model with the lowest AICc and the AICc of the specific model.Bold indicates the most parsimonious model (ΔAICc < 2 and lowest number of parameters).