Nematodes used in the experiments were harvested from monospecific cultures in exponential growth phase, with unidentified bacteria from the habitat as food. Monospecific cultures of four cryptic species of L. marina (Pm I–Pm IV) were raised from a single gravid female per species obtained from the field and maintained on sloppy (1%) nutrient: bacto agar (ratio of 1:4) media (Moens and Vincx, 1998). The nematode stock cultures were incubated in the dark under standardized conditions (temperature: 20°C; salinity: 25) for many generations prior to the start of our experiments. The culture inocula of Pm I, Pm II, and Pm III originate from Paulina salt marsh, Westerschelde Estuary, The Netherlands (51° 21′ N, 3° 49′ E; salinity: 24), whereas Pm IV was obtained from Lake Grevelingen, The Netherlands (51° 44′ N, 3° 57′ E; salinity: 32).

Twenty-five marine bacterial strains (Supplementary Table 1) for the present experiments were obtained from the BCCM/LMG bacterial culture collection at the Department of Biochemistry and Microbiology of Ghent University. All species were selected because they were among the most abundant bacterial taxa found in the microbiome sensu lato of specimens of at least one cryptic species of L. marina in the field (Derycke et al., 2016). Liquid cultures were prepared by isolating single colonies from the stock plates, which were then inoculated in sterile marine broth (Difco) in a 15-ml test tube. The cultures were shake-incubated for 96 h at 180 rpm and 37°C. Prior to this, several test runs and measurements were conducted using a Bio UV visible spectrophotometer at OD₆₀₀ (optical density at 600 nm) to assess the bacterial growth timing. By monitoring the rate of increase in OD_600, we were able to identify the exponential growth phase of the bacterial cultures. Bacterial cell densities were adjusted to 1,000 FAU (Formazin Absorption Units; ISO 7027, 1999), corresponding to an approximate cell density of 10⁹–10¹⁰ cells ml^–1 for each bacterial strain. Considering that OD 600 values are dependent on the size and shape of bacterial cells, and since the bacterial strains used in our experiments differed in cell dimensions, the reported bacterial cell densities represent estimates rather than precise measurements. Bacteria were harvested from the liquid cultures by centrifugation (2,000 rpm; 20-min); the growth medium (supernatant) was removed, and the bacteria (pellet) were washed twice in Artificial Sea Water (ASW; Dietrich and Kalle, 1957) with a salinity of 25 to remove remnants of the growth media and bacterial metabolites, and finally resuspended in ASW. The resulting bacterial suspensions of each species were then diluted in ASW to a density of ∼ 3 × 10⁹ cells ml^–1. This bacterial density supports active population growth and individual fitness of L. marina (dos Santos et al., 2008; De Meester et al., 2011).

Three food diversity treatments (henceforth referred to as “low,” “medium,” “high” diversity) were used in our experiments. Low, medium, and high-diversity mixtures consisted of 5, 10, and 15 bacterial strains, respectively. Replicates of a particular diversity treatment were assigned bacteria by a separate random draw of the appropriate number of species from our pool of 25 bacterial strains. Because we were interested in the effect of diversity, we focused on the number of bacterial strains regardless of their identity. For each treatment, culture aliquots of each component bacteria were mixed in such a way that approximately equal cell numbers of all species were present, and added up to a final density of 3 × 10⁹ cells ml^–1. A suspension of frozen-and-thawed Escherichia coli (strain K12), a known suitable food source for L. marina (Moens and Vincx, 1998; dos Santos et al., 2008; De Meester et al., 2011), was used as a positive control.

A summary of the experimental design for our study is presented in Figure 1.

Parametric tests (ANOVA) were used where possible and non-parametric tests (Kruskal-Wallis for univariate analyses/PERMANOVA for multivariate analyses) were conducted when the assumptions for their parametric alternatives were not met even after transformation. PERMANOVA (Permutation based Multivariate Analysis of Variance; Anderson, 2001) was performed using the Adonis function (Vegan package, Oksanen et al., 2020) based on Euclidean distances with 999 permutations. PERMDISP (Permutational Analysis of Multivariate Dispersions) was always used in conjunction with PERMANOVA to test the homogeneity of multivariate dispersions, thus checking if the observed effect was not the result of variation among the replicates (Anderson, 2004). A post hoc Tukey’s honest significant differences (HSD) test was used for the parametric ANOVA, while pairwise Wilcoxon rank-sum test was conducted for the non-parametric Kruskal-Wallis. Significant factors and interactions for PERMANOVA tests were investigated using posterior pairwise comparisons and p-values were corrected for multiple testing with the Bonferroni method. All statistical analyses were conducted using R (R Core Team, 2021).

Treatment, species, and the interaction between these factors showed significant effects on taxis to food of the four cryptic species of L. marina (all p = 0.001; Figure 3). No significant pairwise differences for the interaction factor remained after Bonferroni correction (all p > 0.05). Higher attraction of nematodes was observed toward all bacterial spots (E. coli, low, medium, high diversity) compared to the negative control (no food) (all p < 0.05). In addition, the nematodes showed higher taxis to medium food diversity compared to E. coli (p = 0.0007), low (p = 0.0007) and high-diversity food treatments (p = 0.01). We also observed significantly higher taxis to food of Pm IV compared to Pm I (p = 0.043) and Pm II (p = 0.001). Within-species analyses further revealed significantly higher attraction of each cryptic species except Pm I toward a medium-diversity bacterial mix compared to the other food treatments (all p < 0.05; Pm I: p = 0.447). In Pm III, high-diversity food also yielded significantly higher attraction than E. coli and low-diversity food (both p < 0.05), and a higher taxis was observed to E. coli as a food source compared to the low-diversity food treatment (p = 0.029).

The cryptic species composition of L. marina in real and fictitious assemblages at every time moment showed very similar patterns for adults and juveniles. Here, we only present the results for adults, which demonstrated that the relative abundances of nematodes were influenced by the interplay between time and interspecific interactions (p = 0.001), and between resource diversity and interspecific interactions (p = 0.002), as well as by the separate effects of time, interspecific interactions and resource diversity (all p < 0.005) (Table 2 and Figure 7).

Both Pm I and Pm III showed higher relative abundances than the other two cryptic species in both fictitious and real assemblages with E. coli as a food source. Pm III also exhibited high relative abundances in fictitious assemblages in low, medium and high-diversity food treatments, yet, it showed the lowest relative abundances among all the cryptic species in real assemblages, i.e., in the presence of interspecific interactions, in all food treatments except the one with E. coli. By contrast, Pm II showed higher relative abundances in the two most diverse food treatments of the (real) multispecies experiment than what could be expected based on monospecific cultures. Pm IV also showed high relative abundances in the presence of interspecific interactions, exhibiting the highest relative abundances of all cryptic species in low and high-diversity treatments. In the medium-diversity food, all cryptic species except Pm III performed well in real multispecies cultures.

Like species composition, the absolute abundances of adults of all nematode species were significantly affected by time, interspecific interactions, resource diversity, and the two-way interactions between these factors (all p < 0.05; Table 3 and Figure 8). Higher absolute nematode abundances were found in fictitious (summed monospecific) than in real assemblages in all food-diversity treatments (all p < 0.05). The lowest abundances were observed with E. coli as the food source, whereas the highest abundances were found in the low-diversity food treatment (all p < 0.05) both in fictitious and real assemblages. The absolute abundances of nematodes were also significantly higher in the high-diversity compared to the medium-diversity food treatment (p = 0.007) in real (multispecies) but not in fictitious assemblages (both p > 0.05).

To our knowledge, the present study is the first to demonstrate taxis to food of nematodes using mixtures of bacterial strains at different levels of diversity. All four cryptic species consistently showed their highest attraction toward the medium-diversity food bacterial mixture, and in three of the four (Pm II, Pm III, and Pm IV), their chemotaxis at this food resource diversity significantly exceeded that toward E. coli and both the low-diversity and high-diversity food. Pm III also exhibited preference toward high-diversity food compared to the E. coli and low-diversity food treatments (Figure 3). These results are particularly striking because different replicates of any given diversity treatments were composed of different combinations of bacterial strains, underscoring the importance of the observed diversity effect. Different bacterial mixtures may elicit distinct bacterial quorum sensing signals and/or produce different types or concentrations of attractants, which may influence the sensory ability of nematodes (Grewal and Wright, 1992; Köthe et al., 2003; Beale et al., 2006). However, in nature, nematodes are exposed to a range of microbiota growing in close proximity to each other and hence to complex mixtures of infochemicals produced by bacteria and other microbial organisms. In this respect, it is relevant to note that the freshwater nematode Bursilla monohystera was attracted to cyanobacterial biofilms and to some of the volatile organic compounds (VOC) they produce. Bursilla monohystera did not respond to any single VOC but a mixture of VOC did elicit a pronounced chemotaxis of the nematode (Höckelmann et al., 2004). Nematodes can sense various signals through several olfactory sensory neurons (Bargmann and Horvitz, 1991; Bargmann et al., 1993), which we hypothesized may function optimally under conditions where a sufficient diversity of stimuli is present, but suboptimally when too large a mixture of slightly different stimuli has to be “processed”. If so, this might explain why cryptic L. marina species were more attracted to mixtures of 10 bacterial strains than to lower and higher-diversity mixtures. Another interesting aspect which remains largely unexplored is the relationship between marine nematodes and fungi, and the potential of the latter as food sources and/or as attractants for nematodes. Although fungal-feeding nematodes are lacking in marine nematode benthic assemblages, except for some algae-associated Aphelenchoididae (Moens et al., 2013), there is evidence from barcoding studies of some association between fungi and marine nematodes (Bhadury et al., 2009, 2011; Bhadury and Austen, 2010) indicating that a fungal microbiome indeed exists in these organisms.

But did the observed taxis-to-food patterns concur with any fitness responses of these nematodes to resource diversity? We evaluated both population-level and individual nematode fitness, the latter on the basis of two traits which both contribute to the rate of population increase, i.e., minimum juvenile development time and fecundity. In general, there was a trend of shorter juvenile development time (Figure 4) and higher fecundity (Figure 5) at high and medium-diversity food; however, only some of these effects were statistically significant (Table 4). In addition, we observed a higher fecundity of Pm III compared to the other cryptic species, in accordance to De Meester et al. (2015a), yet, our study revealed that this was more pronounced in the high-diversity food treatment.

At the population level, each bacterial mixture supported population growth of the nematodes (preliminary analyses also showed that the individual bacterial strains can support nematode growth), but a more diverse food source enhanced the population growth of L. marina, especially for Pm I and Pm IV (Table 4 and Figure 6). Pm II and Pm III exhibited a similar response, but performed better in all bacterial mixtures (low, medium and high-diversity food) than on the single-strain resource E. coli. Our findings appear to corroborate the “balanced diet” hypothesis, which suggests that a more diverse resource offers a more complete and/or complementary range of nutritional resources, leading to a higher production of consumer biomass (DeMott, 1998). By contrast, both the “variance-in-edibility” hypothesis and the “dilution hypothesis” or “resource concentration hypothesis” predict a negative impact on consumers when resource diversity is increased and are supported by a considerable body of evidence (e.g., Andow, 1991; Steiner, 2001; Wilsey and Polley, 2002; Hillebrand and Cardinale, 2004; Duffy et al., 2005), whereas examples of the balanced diet effect of diversity on consumers are more scarce. Nevertheless, several studies have demonstrated enhanced herbivore growth and biomass accumulation with mixed diets of primary producers than with single-species diets (e.g., Pfisterer et al., 2003; Worm et al., 2006). However, in the rapidly changing and ephemeral habitat where L. marina lives, the coexistence of the cryptic species is likely not only determined by resource availability, but also by the responses of the cryptic species and their microbiomes to environmental changes (Derycke et al., 2006; De Meester et al., 2015a) and/or their ability to disperse to suitable patches with more favorable conditions (De Meester et al., 2012), such as a site with high resource diversity.

Resource diversity can also affect the response of cryptic species to competition. We found that competitive interactions between the cryptic species of L. marina were substantial, as evidenced by the fact that when all four cryptic species were inoculated together, the total abundance of nematodes was lower than the sum of the abundances of the four species in the monospecific experiments (Figure 8). These competitive interaction effects corroborate the findings of previous research (De Meester et al., 2011, 2015b). In addition, the total numbers of nematodes were lower with E. coli as food compared to a more diverse food source both in real and fictitious assemblages, showing that nematode abundances of the mixed species assemblages were also clearly influenced by resource diversity. Based on the principle of competitive exclusion, species occupying the same ecological niche cannot stably coexist due to strong interspecific competition for resources (Volterra, 1926; Gause, 1936). Niche diversification, for instance through resource partitioning, may alleviate competition (Sánchez-Hernández et al., 2017).

Interspecific interactions also had an effect on the species composition of the nematode assemblages. In real assemblages (multi-species experiments), Pm I was the competitively strongest species in the E. coli food treatment in which Pm IV almost became excluded after 21 days. Pm I and Pm IV are the two most closely related cryptic species and they seldom co-occur in the field (Derycke et al., 2006, 2008), so we expected substantial competition between them. The present results also corroborate that when simultaneously inoculated in a homogeneous environment without spatial structure and with E. coli as a food source, Pm I is the strongest competitor, followed by Pm III, while Pm II and Pm IV are poor competitors (Figure 7).

However, the competitive hierarchy among these four species is not fixed. All four species can co-exist for longer time periods when dispersal opportunities are present (spatial structure) (De Meester et al., 2014). In addition, the outcome of the competition depends on the precise species combination and on environmental factors such as salinity and temperature (De Meester et al., 2011, 2015a). The present study demonstrates that resource diversity also affects the composition of assemblages of cryptic L. marina species. In strong contrast with the E. coli treatment, Pm IV exhibited the highest relative abundances of all species in low and high-diversity food treatments in multispecies experiments; Pm I and Pm II also occurred in high abundances. On the other hand, Pm III exhibited the lowest relative abundances at all resource diversity levels except the E. coli treatment in multispecies experiments (real assemblage), although it exhibited the highest relative abundances in fictitious assemblages, demonstrating that this species suffered from interspecific interactions irrespective of resource diversity (Figure 7). Derycke et al. (2016) argued that Pm I and Pm III partly feed on different bacterial strains, and that Pm III is a much less generalist feeder than Pm I. Since monospecific cultures of Pm III performed well on each of the bacterial diets offered here, suitable food for Pm III was sufficiently available in all food regimes. Given the poor performance of Pm III in the low, medium and high-diversity food treatments, other species were likely more successful at utilizing these resources than Pm III. When a competitively dominant species expands its niche, the niche space available for other species is reduced, resulting in selection pressure to the inferior competitor to use whatever resources are left available (Levins, 1968; Jasmin and Kassen, 2007; Pekkonen et al., 2013). In this case, there may be a pronounced effect of competition on selective feeders since their preferred food may become limiting in the presence of other cryptic species, whereas more generalist feeders may be able to minimize competition in a more diverse food environment through enhanced resource partitioning. When resource diversity is high, the competing species may segregate in resource utilization by specialization, which may minimize the negative effects of competition (Sánchez-Hernández et al., 2017).

In the field, the local diversity of marine nematodes is typically high (Heip et al., 1985; Moens et al., 2013). Despite their diversity, quantitative and qualitative information on the feeding ecology of nematodes at the species, genus and even family level, remain scant (Montagna, 1984; Moens et al., 2014; Wu et al., 2019). Marine nematodes have traditionally been classified by the morphology of their stoma, and, to a lesser extent, based on pharyngeal musculature (Wieser, 1953). However, our findings show that even very closely related nematodes with highly similar stoma, such as the cryptic species of L. marina, can have differential responses to resource diversity, which can impact the outcome of interspecific interactions. This indicates that the link between morphology and feeding, although informative, is neither as absolute nor as straightforward as proposed by Wieser’s classification (Moens et al., 2004), which has significant implications from a functional perspective and in understanding how resource partitioning affects the composition and diversity of nematode communities.