The P. hermaphrodita strain used in this study (DL 309) was isolated from a dead slug collected in Oregon, USA (Howe et al., 2020). Nematodes were maintained in the lab on standard nutrient growth media (NGM) agar plates with co-cultured bacteria as a food source. NGM plates were prepared using autoclave-sterilized media, equipment, and sterile plastic plates (VWR International, Lutterworth, UK) and seeded with the appropriate bacterial cultures. Prior to our study, nematodes were bleached with sodium hypochlorite 0.25 M solution following standard worm bleaching protocol (Stiernagle, 2006) to minimize bacterial carry-over from previous culturing conditions. Nematodes were monitored on plain NGM agar plates for 24 h to ensure no bacterial growth, then harvested and washed once with sterile M9 buffer and twice with molecular biology grade water. Next, nematodes were transferred onto fresh NGM plates (Stiernagle, 2006) seeded with one of the three designated bacterial food sources: E. coli OP50 (EC), Pseudomonas sp. (PS)., or the original bacterial community (BC) associated with the nematodes. The E. coli OP50 culture was obtained from the Caenorhabditis Genetics Center at the University of Minnesota. The Pseudomonas strain was isolated from the original bacterial community that co-cultured with a P. hermaphrodita sampled at the Oregon State University campus (Mc Donnell et al., 2018). We confirmed the genus identity of the strain with 16S rRNA and rpoB gene sequencing.

Prior to DNA isolation, nematode samples were homogenized thoroughly by bead beating. We extracted DNA from the homogenate and from the negative control samples using the PowerSoil DNA Extraction kit (MOBIO, Carlsbad, CA, USA). DNA extractions were then quantified using a fluorescent plate reader at the OSU Center for Genome Research and Biocomputing (CGRB).

16S library preparation, amplification, and sequencing were conducted by the OSU CGRB. Amplicon libraries were prepared following the standard protocol by Illumina (Illumina, Inc., 2013). Primers targeting the V3-V4 region (forward: 5’TCGTCGGCAGCGTCAGATGTG TATAAGAGACAGCCTACGGGNGGC WGCAG3’; reverse: 5’GTCTCGTGGGCTCGGAGATGTG TATAAGAGACAGGACTACHVGGGT ATCTAATCC 3′) were used for PCR amplification (Klindworth et al., 2013). Subsequently, the PCR products were processed through clean-up, index PCR, library quantification, normalization, and pooling. Libraries were sequenced using paired-end MiSeq v2 Nano 300 bp platform (Illumina, USA).

We used the DADA2 pipeline (Callahan et al., 2016) to trim reads, merge paired reads, denoise, quality filter, and infer amplicon sequence variants (ASVs). ASVs represent exact 16S rRNA gene sequence variants resolved to the level of single-nucleotide dissimilarity. Compared to the commonly used operational taxonomic units (OTUs) with a threshold of nucleotide difference of 97%, ASVs have been shown to offer better reproducibility, reusability, and comprehensiveness in 16S microbiome analyses (Callahan et al., 2017).

Phylogenetic analyses and representations of the ASVs were constructed using MAFFT alignment (Katoh et al., 2002) and FastTree 2 (Price et al., 2010). To assign taxonomy to the ASVs, we utilized a self-trained classifier, trained on the Silva Project’s database (release 132) (Quast et al., 2013). Putative contaminant sequences were identified and filtered out of sample data using the R package decontam (Davis et al., 2018). ASVs were considered a contaminant if they were more prevalent in negative controls than in positive samples and/or their frequency significantly varied inversely with sample DNA concentration (probability threshold p < 0.1).

We utilized the Phylum, Genus, and ASV levels as units for subsequent taxa analyses. Since different taxonomic classifiers trained on different databases may differ at lower taxonomic ranks (Balvočiūtė and Huson, 2017), the Phylum level was chosen as the main unit for taxa analysis; however, the Genus and ASV level were also selected for a finer resolution and subsequent composition analysis.

We used four metrics: Chao1, Shannon Index, Inverse Simpson, and Faith’s Phylogenetic Diversity (PD) to compare and evaluate α-diversity from different approaches. Chao1 index is an estimator of species richness, i.e., the expected total number of OTUs/ASVs in the sample given all the species were identified. It weighs the low abundance taxa (i.e., only singletons and doubletons) to infer the number of missing species (Chao, 1984). The Shannon Index estimates the overall richness and also the evenness/uniformity between the taxa present in the sample (Spellerberg and Fedor, 2003). Also focusing on the species richness and evenness, however, the Inverse Simpson index emphasizes more on species evenness (Simpson, 1949), while the Shannon Index gives greater weight to the species richness. Faith’s PD metric considers the phylogeny of taxa to estimate diversity and is proportional to how much of the OTU/ASV phylogenetic tree is covered by the taxa present (Faith, 1992).

To calculate the Shannon Index and Inverse Simpson, the ASV table was rarefied ten times to a depth of 100,000 reads to lessen the random biases introduced by random subsampling. Rarefying the ASV table to a fixed number of reads per sample is necessary for calculations of the Shannon Index and Inverse Simpson, as these two abundance-based indices could be substantially affected by the differences in the total number of reads between samples. The two indices were calculated for each of the ten rarefied tables using the package phyloseq (McMurdie and Holmes, 2013), and the mean values were used for α-diversity analysis. As for Chao1, which depends on low-abundance taxa, the ASV table was not rarefied to retain all ‘rare’ taxa. The Faith’s PD index was calculated using the package picante (Kembel et al., 2010). We assess the significance of the differences among the bacteria-nematode combinations at the same time point using Kruskal-Wallis tests. Mann–Whitney U tests or paired t-tests were applied, as appropriate, to test for the difference between the two time points before and after infection. The normality of diversity values was assessed using the Shapiro–Wilk normality test.

To assess the difference in the microbial community between samples, we used four β-diversity metrics: two abundance-weighted β-diversity indices (Bray–Curtis dissimilarity and weighted UniFrac) and two indices based on the species (ASV) presence-absence matrix (Jaccard distance and unweighted UniFrac). While Bray-Curtis and Jaccard are non-phylogenetic metrics, weighted and unweighted UniFrac also take into account the phylogeny of the bacterial community composition to assess differences between samples. Specifically, UniFrac distances are based on the sum of branch length shared between samples on the phylogenetic tree constructed from all the 16S rRNA sequences from all communities of interest. To test whether the groups differ significantly from each other, we conducted the permutational multivariate analysis of variance (PERMANOVA) test using the package vegan (Oksanen et al., 2012).

Slugs that were exposed to the PS high dose treatment (LT50: 5.68 days) survived for significantly (p < 0.05) fewer days than those treated with both the EC high (LT50: 7.74 days) and low (LT50: 8.61 days) rates (Figure 1). For BC high (LT50: 5.92 days), slugs survived significantly (p < 0.05) fewer days than BC low (LT50: 6.77 days) and both EC rates. Lastly, slugs exposed to BC low survived significantly fewer days than those treated with EC high (Figure 1).

After quality filtering and removing chimeric sequences, we obtained a total of 10,485,882 reads, with an average of 582,549 reads and a median depth of 349,960 reads per sample. 2,938 ASVs (i.e., taxa) were identified. The decontamination procedure removed 131 ASVs, leaving 2,807 unique ASVs for subsequent analyses. The components of the removed contaminants are shown in Supplementary Figure S1. The filtered data set then contained 10,422,195 reads in total, with an average depth of 579,010 reads and a median 343,942 reads per sample. Details of the number of quality-controlled reads obtained in each sample, as well as the total number of reads, median, and mean values of reads per nematode-bacterial combination were listed in Supplementary Table S2.

The number of different genera detected in each sample included in the analysis is listed in Table 1. On the Phylum level, the majority of bacterial components in all samples were identified as members of the phyla Bacteroidetes and Proteobacteria. We observed a shift in composition towards Proteobacteria in the BC-PostInf and PS-PostInf samples compared to their PreInf counterparts, while in the EC-PostInf samples, there was a shift towards Bacteroidetes (Figure 2A). In general, the number of genera in all samples increased almost fourfold after infection. The abundance of the main genera also altered remarkably between PreInf and PostInf samples (Figure 2B).

In P. hermaphrodita’s original microbiome (i.e., the BC-PreInf community), 15 phyla were identified with Bacteroidetes being the most abundant phylum, followed by Proteobacteria. On the Genus level, the most abundant groups were Pseudochrobactrum, Flavobacterium, Raoultella, and Pseudomonas out of the 82 genera observed (Figure 2B). Details of relative taxa abundance on the Phylum and Genus level in the BC-PreInf communities were shown in Supplementary Table S3A.

The numbers of ASVs unique to specific sample sets and shared among sample sets were illustrated in Figure 3. In general, the number of ASVs unique to a single sample set was higher than those shared between different sets. The BC-PostInf samples had the largest number of unique ASVs (i.e., 600) that were not shared with any other type of samples, followed by the EC-PostInf and PS-PostInf (451 and 337 ASVs, respectively). Comparing PreInf and PostInf samples, the BC-PreInf and BC-PostInf communities shared 82 common ASVs, while PS-PreInf and PS-PostInf shared 36 ASVs, and 44 ASVs were shared between EC-PreInf and EC-PostInf samples.

We hypothesized that slug infection would alter the diversity of the bacterial community within P. hermaphrodita nematodes, and that these changes would be influenced by the pre-infection bacterial associate. To test this, we evaluated α-diversity of the nematodes’ bacterial communities using four metrics: an estimator of species richness (Chao1), two metrics of species richness and evenness (Shannon index, Inverse Simpson), and a phylogeny-based measure (Faith’s Phylogenetic Diversity - PD) (Figure 4).

Bacterial community diversity significantly increased following nematode infection when comparing PreInf and PostInf samples across all three bacterial combinations. Specifically, the Chao1 index showed a significantly higher median across all PostInf samples compared to PreInf (Mann–Whitney U test, p < 0.001) (Figure 4A), indicating a clear increase in the number of ASVs and higher total species richness within the PostInf samples. Faith’s PD also showed trends toward higher diversity in terms of species richness and phylogenetic diversity post-infection, though this difference was only suggestive (Mann–Whitney U, p = 0.06).

Examining each treatment individually, bacteria-nematode combinations exhibited distinct diversity patterns following infection (Figure 4B). The EC microbiomes showed significant increases in Chao1 richness, Inverse Simpson index, and Shannon diversity following slug infection (Wilcoxon test, p = 0.05). These changes were consistent with each other and indicated a higher total number of bacterial taxa (species richness) and higher evenness among taxa members post-infection, though insignificant Faith’s PD result suggested that there was no evidence of a shift in phylogenetic diversity of taxa compositions. In BC PostInf samples, Faith’s PD and Chao1 richness increased significantly (Wilcoxon test, p = 0.05), while the Inverse Simpson index decreased (Wilcoxon test, p = 0.05), demonstrating an increase in phylogenetically diverse taxa alongside reduced evenness, which possibly resulted from the emergence of many new low-abundance taxa compared to pre-infection samples. For the Pseudomonas-enriched communities, the PS-PostInf samples demonstrated strong evidence of increases in Shannon diversity and Chao1 index (Wilcoxon test, p = 0.05) compared to PS-PreInf, indicating that the bacterial communities in the Pseudomonas-enriched treatment developed greater species richness and community diversity following slug infection.

As expected, α-diversity differed significantly among microbiomes from the three bacteria-nematode combinations prior to infection (Kruskal-Wallis test, p = 0.05), with EC-PreInf and PS-PreInf communities showing lower diversity than BC-PreInf due to their more defined, single taxon-enriched starting compositions. However, after the infection assay, α-diversity indices between the three bacteria-nematode combinations did not differ significantly (Kruskal-Wallis test, p = 0.30). This suggests that the infection process introduced additional taxa to the bacterial communities, potentially reaching a similar α-diversity equilibrium regardless of the initial bacterial associate.

Given the distinct starting bacterial associates, we hypothesized that bacterial community structures would differ significantly among the three bacteria-nematode combinations (EC, PS, and BC). We also predicted that going through the slug infection process would alter the bacterial community composition within each treatment, leading to detectable shifts in community structure between pre- and post-infection samples. To evaluate the dissimilarity in bacterial community composition among samples from the three combinations before and after the infection assay, we applied four complementary β-diversity indices: an abundance-weighted metric (Bray–Curtis dissimilarity), a presence-absence index (Jaccard distance), and two phylogenetically-based indices with and without weighting on taxa abundances (weighted and unweighted UniFrac, respectively).

Bray-Curtis analysis results showed a clear pattern of distinct sample clustering (Figure 5). The bacterial community structures were significantly different among the EC, PS, and BC sample sets (PERMANOVA, p = 0.001, R2 = 0.541). This result was robust to Jaccard distance (PERMANOVA, p = 0.001) and weighted UniFrac indices (PERMANOVA, p = 0.001). However, unweighted UniFrac showed only marginal significance (PERMANOVA, p = 0.06). These results suggest that microbiomes from distinct bacterial conditions significantly differed in species composition, phylogenetic distribution, and community structure. Since unweighted UniFrac considers only the presence or absence of taxa and is particularly sensitive to differences in rare lineages, while weighted UniFrac incorporates taxa abundance information, the significant differences shown by weighted UniFrac but not unweighted UniFrac suggest that the differences among the three bacterial treatments were mainly driven by variations in the relative abundances of common taxa, rather than differences in the presence or absence of rare, phylogenetically unique species.

Comparing pre-infection and post-infection samples, Bray–Curtis dissimilarity, Jaccard distance, and unweighted UniFrac indicated a significant compositional difference (PERMANOVA, p = 0.02, p = 0.01, p = 0.001, respectively); however, weighted UniFrac did not detect significant changes (p = 0.25). The shift detected by unweighted UniFrac but not by abundance-weighted UniFrac suggests that the differences in community composition and structure pre- and post- infection mainly resulted from changes in rarer taxa, which contribute less to the abundance-weighted phylogenetic distance calculations.

We analyzed the ASV components of the PS microbiome, with specific focus on those deriving from Pseudomonas spp., to evaluate the potential effects of the infection process on this group of bacteria (Figure 6). The microbiome of PS-PostInf samples only shared a total of 36 ASVs with PS-PreInf, four of which were identified as belonging to the genus Pseudomonas. The relative abundance of these four ASVs in PS-PreInf samples was 13.23% of the total, and expanded 4.5-fold to 59.58% in PS-PostInf samples. Among these four ASVs, one was present in all three bacterial-nematode combinations sampled both before and after infection, and displayed an increase of 7.5, 3.6, and 3.2% in EC, PS, and BC communities, respectively. The three other Pseudomonas ASVs were detected only in the microbiomes of PS and BC nematodes, both before and after the infection assay; they were absent from all EC samples. While their abundances all increased after the course of infection in PS samples - by 6.9, 30.5, and 5.5% - there was a consistent decrease of these ASVs by 1.38, 4.13, and 0.73% from BC-PreInf to BC-PostInf samples. Two Pseudomonas ASVs present exclusively in PS-PreInf accounted for just 0.004% abundance; fourteen ASVs unique to the PS-PostInf samples had a total relative abundance of 1.11%.

We also hypothesized that nematode-associated microbiomes would change significantly from pre-to post-infection samples. Host infection might act as a selection pressure on the bacterial community, favoring the bacteria that can better contribute to the nematode’s virulence, spread, reproduction, and killing efficacy. Furthermore, recent research has shown differential microbiome responses between nematode-susceptible and nematode-resistant host gastropod species (Sheehy et al., 2022).

The microbiome samples of all the EC, PS, and BC samples showed complex taxa (ASVs) compositions. The species richness and diversity of EC-PreInf and PS-PreInf communities were higher than our expectations, given that the nematodes had been bleached with hypochlorite, washed thoroughly with water and salt solutions, and reared on pure cultures of E. coli or Pseudomonas sp. before the infection assay. The intended bacteria were present in the samples collected before infection trials; however, Pseudomonas sp. and E. coli accounted for a relatively minor component of the community by the time the infection experiments began. Several factors likely contributed to the unexpectedly high taxonomic diversity observed in the EC and PS samples, including environmental bacteria in the lab, certain bacteria surviving the bleaching procedures, and the possibility that there are unknown bacterial endosymbionts associated with P. hermaphrodita as has been reported in other nematode systems (Brown et al., 2018; Mobasseri et al., 2019).

Before the infection trial, fifteen different phyla and 82 genera were observed in the BC-PreInf samples, with Pseudochrobactrum, Flavobacterium, Raoultella, and Pseudomonas being the four most abundant taxa in the community. These four genera have all been reported in association or rivalry with nematodes, for example, Pseudochrobactrum spp. with the EPN Steinernema (Ogier et al., 2020), Raoultella spp. with the root-knot and pine wood nematodes (Liu et al., 2019), Flavobacterium spp. with the insect-parasitic nematode Rhabditis blumi (Park et al., 2011), the soybean cyst nematode Heterodera glycines (Nour et al., 2003), the oriental beetle pathogenic nematode Butlerius spp. (Yi et al., 2007), and Pseudomonas spp. with C. elegans (Pees et al., 2024). The high-abundance components in BC-PreInf samples are relatively similar to those of C. elegans, whose dominant genera are Pseudomonas, Stenotrophomonas, Ochrobactrum, Sphingomonas, and unclassified Enterobacteriaceae (Raoultella is a representative of this family) (Dirksen et al., 2016; Zhang et al., 2017). Stenotrophomonas, Ochrobactrum, and Sphingomonas spp. were all detected in the top 10 most abundant genera in our BC-PreInf bacterial community.

The complex nature of the bacterial community structure and the results of slug mortality during the infection assay indicated that P. hermaphrodita could carry out slug infection with various bacterial partners regardless of the primary bacterial food source, although at different rates. This observation agrees with a previous finding that the nematode is able to associate with non-specific and complex microbial assemblages (Rae et al., 2010). However, slug-killing efficacy clearly varied between different bacteria-nematode combinations. BC-cultured nematodes were the first out of the three treatments to cause significant slug mortality on Day 5, followed by the PS nematodes on Day 6, and lastly, the EC nematodes, starting Day 10. However, the nematodes treated with Pseudomonas took only 3 days to kill all slugs and were the first to cause 100% slug mortality, followed by BC nematodes (5 days) and EC nematodes (10 days). This accelerated killing rate suggests that Pseudomonas spp. may influence the nematode’s pathogenicity, possibly by assisting P. hermaphrodita’s virulence, acting as an opportunistic pathogen, and/or providing a favorable food source that benefits the nematode’s growth and development, thereby enhancing infection and killing efficacy.

Analysis of ASVs in both the PS-PreInf and PS-PostInf sample sets from the Pseudomonas treatment revealed four common ASVs belonging to the genus Pseudomonas (see Results). In the PS treatment, the relative abundance of these four Pseudomonas ASVs expanded 4.5-fold from PS-Pre-Inf to PS-PostInf samples, accounting for approximately 60% of total abundance of all taxa in the community. Among these four ASVs, we detected the presence of one particular ASV in all sample sets, whose relative abundance consistently increased after the course of infection in all three bacterial treatments. This suggests that this specific ASV might be stably retained by the nematode and could be important to the slug infection activity. The abundance of the three other Pseudomonas ASVs also expanded after infection in PS-PostInf samples. However, in BC-PostInf samples, Pseudomonas ASVs did not exhibit the same expansion pattern. The lack of Pseudomonas expansion in BC-PostInf samples, despite the apparent pathogenic effects, suggests that the potential role of Pseudomonas spp. in pathogenicity may be dependent upon the microbiome context. In the PS treatments, which have lower species complexity than BC, the expansion of these ASVs may suggest a role in pathogenicity. This pattern is not observed in BC, potentially due to interactions or competition with other bacterial taxa that interfere with Pseudomonas function in the infection process, preventing it from achieving the dominance necessary to enhance nematode pathogenicity as observed in the Pseudomonas-enriched cultures.

In general, we found that α-diversity indices were higher in all post-infection samples compared to pre-infection samples, especially in the Pseudomonas and bacterial community treatments. For BC-PostInf samples, Chao1 and Faith’s PD α-diversity indices increased significantly, which may indicate the appearance of low-abundance taxa during the course of infection, and also suggest that the diversity in the BC PostInf microbiome results not just from the presence of a few highly diverse taxonomic groups, but rather from a broader range of phylogenetically different taxa. Indeed, for example, the total number of different genera in the BC microbiomes more than tripled after infection, however, the composition of the top 20 most abundant taxa remained stable. Likewise, for the Pseudomonas-enriched treatment, we reported significant increases in Chao1 and Shannon indices after the assay, indicating a spike in rare taxa post-infection, although composition analysis pointed to the dominance of Pseudomonas taxa in the community.

We found that the difference between samples from different bacteria-nematode combinations was significant for Bray–Curtis dissimilarity, Jaccard distance, and weighted UniFrac metrics, but not for unweighted UniFrac. This reveals that between bacterial combinations, the communities comprehensively differed in structure, taxa composition, and phylogenetic distribution. The non-significant unweighted UniFrac result suggests that the same major phylogenetic groups were present across all bacterial treatments. However, the significant weighted UniFrac result implies that the relative abundance of these shared phylogenetic lineages was drastically different. In other words, the microbiomes in different bacterial combinations were not reshaped by the loss or gain of entire phylogenetic branches, but rather by profoundly altering which taxa were dominant.

Meanwhile, samples collected before and after infection show a significant difference according to the unweighted UniFrac but not the weighted UniFrac metric. This pattern indicates that the shift in community composition pre- and post-infection mainly resulted from the gain or loss of low-abundance taxa, while the most abundant taxa remained relatively stable. This observation is in concordance with the findings of the α-diversity analyses discussed above.