Microbiome Differentiation Among Coexisting Nematode Species in Estuarine Microhabitats: A Metagenetic Analysis

1 Introduction

Coexistence of functionally similar and/or taxonomically related species has long intrigued ecologists, because species with almost identical ecological niches/functions cannot coexist under a purely deterministic coexistence model since their intense competition would ultimately lead to exclusion (Hardin, 1960). Under the ‘niche concept’, natural selection drives the evolution of species into divergent and complementary niches (Hector and Hooper, 2002; Chase and Leibold, 2003). The more divergent their niches, the less intensely species will compete, and the more prominently intraspecific rather than interspecific competition regulates population dynamics (Leibold and McPeek, 2006). Niche differentiation insights are based on the differential tolerances of species to abiotic conditions and/or to their differential resource utilization, which can affect demographic traits and population fitness, resulting in a higher competence of one species over others, given a particular niche (Chase and Leibold, 2003). A stable species coexistence can therefore be maintained when co-occurring species are favored under different sets of environmental conditions, inhabit different structures of the same habitat (i.e., microhabitats), and/or systematically utilize different resources (Abrams, 1983; Ekschmitt and Griffiths, 1998; Chesson, 2000). None of the above is absolute, particularly when small spatial scales are considered, as species share very similar conditions, microhabitats and resources in the same habitat at the same time. Despite such niche overlap, differentiation of the ‘realized niche’ of closely related species under competition can occur, leaving some ‘space’ within the niche for competent species to persist (Hutchinson, 1959; Bruno et al., 2003). Conversely, when species with (almost) identical ecological traits and fitness compete locally (Chesson, 2000), their coexistence can be attained due to neutral dynamics (Hubbell, 2001).

In the sediments of aquatic ecosystems, nematodes are often the most abundant and diverse invertebrates (Heip et al., 1985; Traunspurger, 2000; Moens et al., 2013). A large number of species can co-occur at very small spatial scales (e.g., at submeter scales), holding a number of horizontal interactions, including competitive, facilitative and inhibitory responses (De Mesel et al., 2006; dos Santos et al., 2009; dos Santos and Moens, 2011; De Meester et al., 2015a). Based on commonly used functional guild classifications (Wieser, 1953; Bongers et al., 1991), closely related species (confamiliar, congeneric) are implicitly assumed to be confunctional. This has been recently overturned by the discovery of cryptic species within the same morphospecies complex that have sympatric distributions (Derycke et al., 2005; Derycke et al., 2007a; Derycke et al., 2012; Derycke et al., 2013; Ristau et al., 2013) and non-identical functions (De Meester et al., 2016) and exhibit differential tolerances to environmental conditions (De Meester et al., 2011; Van Campenhout et al., 2014; De Meester et al., 2015a; De Meester et al., 2015b; Monteiro et al., 2018). Evidence for microhabitat and resource differentiation among cryptic nematode species has also been provided (Derycke et al., 2016; Guden et al., 2018). Nevertheless, knowledge on niche divergence and its mechanisms between sympatric nematode species belonging to the same functional guild is scant (Venette and Ferris, 1998; Blanc et al., 2006; Estifanos et al., 2013), particularly in marine nematodes. For example, bacterivorous nematodes of the families Monhysteridae and Rhabditidae, each usually represented by several species, frequently co-occur within various microhabitats (e.g., decomposing macroalgae) in the intertidal zone of coasts and estuaries (Derycke et al., 2006; Derycke et al., 2007a; Derycke et al., 2008; dos Santos et al., 2009). At least some degree of resource differentiation has been suggested among confamiliar species of Monhysteridae (Moens et al., 1999) and among cryptic species of Litoditis marina (Derycke et al., 2016) originating from the same location and microhabitats, but the extent of such resource differentiation remains poorly understood. Furthermore, the decomposition status of macroalgae has been suggested to play a role in substrate colonization preferences of closely related nematode species (Moens et al., 1999), which may reflect resource differentiation in relation to decomposition.

High-throughput sequencing technologies (HTS) nowadays allow the characterization of microbial communities, including remnants of ingested bacteria, associated with bacterial-feeding nematodes (Derycke et al., 2016; Schuelke et al., 2018; Rzeznik-Orignac et al., 2018). Microbiomes normally comprise the genetic diversity of all microorganisms associated with a host, including the gut microflora (commensal bacteria), other (endo-/ecto-) symbiotic or pathogenic bacteria (called microbiome ‘sensu stricto’), and bacteria utilized as food or attached to it (microbiome ‘sensu lato’) (Musat et al., 2007; Derycke et al., 2016; Dirksen et al., 2016; Schulenburg and Félix, 2017). Therefore, microbiome characterization has potential to explain niche partitioning linked to resource availability and utilization, but also linked to symbiotic relationships with microbial organisms (Schuelke et al., 2018).

Here we aim to explore differences in microbiomes of bacterivorous nematode species that co-occur in the same habitats to provide insights into resource use and niche differentiation of functionally similar species. Co-occurring monhysterid and rhabditid nematodes from microhabitats in the intertidal zone of the Scheldt Estuary, the Netherlands, were collected for individual microbiome characterization, using amplicon sequencing of the bacterial 16S rRNA gene. Based on the classical coexistence theory and the niche differentiation concept, we expected that resource partitioning would be more prominent among species, thus resulting in greater microbiome differences within shared microhabitats and less among microhabitats. Furthermore, we anticipated that nematode species would exhibit some degree of plasticity in their resource use and thus, they would have partly different microbiomes in different microhabitats, a main factor for niche divergence that could contribute to microhabitat diversification. We also hypothesized that resource differentiation of species within microhabitats may vary in time, along with the natural temporal variability of the microbial community associated with algae and cordgrass in the field, or with decomposition status of the phytal microhabitat where the nematodes occur.

2 Materials and Methods

2.1 Sampling

Nematodes were sampled by collecting fractions of their potential microhabitats (here macroalgae and saltmarsh plants) in the Paulina salt marsh in the polyhaline reach of the Scheldt Estuary, the Netherlands (51°21’01.9”N 3°43’09.1”E). Samples of algae or saltmarsh material were collected at low tide at two sampling occasions representing two sampling seasons, in November 2015 (autumn) and in April 2016 (spring). Microhabitats were chosen as a function of the ‘target’ nematode species (Table 1) we were aiming to collect, based on information from a year-round monthly survey of the co-occurrence of abundant bacterivorous species in this saltmarsh (Moens et al., 1999). The sampled microhabitats included: decomposing thalli of the brown algae Fucus vesiculosus L. and Fucus spiralis L., fresh and decomposing thalli of the green alga Ulva sp. L., and decomposing but still attached as well as detached leaves of the cordgrass Spartina anglica C.E. Hubb. Decomposing Ulva thalli and old/attached Spartina leaves were only found in November. Three replicates per type of material (i.e., an entire algal thallus or, in the case of fragmented material or cordgrass leaves, an approximately equal amount per replicate) were haphazardly collected from the same location (within 1m scale for microhabitats of the same type and within few meters at the same site for different microhabitat types, e.g. cordgrass and macroalgae) on each sampling occasion, to increase the chance of obtaining enough specimens of the ‘target nematode species’. These were transferred on ice and stored frozen in the laboratory at -20°C until further processing (< 2 months after sampling). Samples of algae or cordgrass were washed carefully into a 5L glass beaker with distilled water, which was then decanted and sieved over a 32-μm mesh, and the supernatant was kept in smaller vials for direct observation of the nematodes. This procedure was repeated 5 times for each sample.

Table 1 Identified nematode species and number of specimens collected from each microhabitat type in autumn and spring.

2.2 DNA Extraction and Nematode Identification

Morphology-based identification and DNA barcoding was used to create 345 distinct samples of nematode specimens collected from 30 algal/plant material samples. Nematodes were first identified under an inverted microscope into 10 species belonging to 3 families (i.e., Rhabditidae, Panagrolaimidae, Monhysteridae) (Table 1), and individuals of the ‘target species’ were sorted in a counting dish and hand-picked with a fine Tungsten wire under a stereomicroscope for identification to species level using DNA-barcoding. They were bathed one by one in an embryo dish with molecular-grade water to remove adhering particles and bacteria. Each specimen was transferred to 20 μl of Worm Lysis Buffer (WLB) (50 mM KCl, 10 mM Tris pH 8.3, 2.5 mM MgCl₂, 0.45% NP40 (Tergitol Sigma), 0.45% Tween20) in a single microcentrifuge tube (Eppendorf AG, Hamburg, Germany), and stored in the freezer (at -20°C) before DNA extraction. DNA was extracted by adding 1 μl of proteinase K (10 mg ml^-1) in the tube, which was incubated for 60 min at 65°C, then for 10 min at 95°C, and finally centrifuged for 1 min at maximum speed (13 000 rpm). DNA extracts were stored at -20°C until further processing.

1 μl of the extracted DNA per sample was used as template DNA for 24 µl of PCR mixture (containing 2.5 µl of 10 × PCR buffer (15 mM MgCl₂ (Qiagen), 2.5 µl loading dye (Qiagen), 2 µl MgCl₂ (25 mM), 0.5 µl dNTP (10 mM), 0.125 µl of each primer (25 µM), 0.125 µl TopTaq DNA polymerase (Qiagen), and 16.125 µl PCR grade water) for PCR amplifications. The primer set G18S4 forward (5’-GCTTGTCTCAAAGATTAAGCC-3’) (Blaxter et al., 1998) and 4R reverse (5’-GTATCTGATCGCCKTCGAWC-3’) (Armenteros et al., 2014) was used to amplify part (ca. 800 bp) of the 18S rRNA gene (small subunit or SSU). PCR amplification was performed with the following optimized cycle conditions: an initial denaturation step for 5 min at 94°C, 40 cycles of denaturation for 30 s at 94°C, annealing for 30 s at 56°C, and extension for 1 min at 72°C, and a final extension step for 10 min at 72°C. PCR products were loaded on agarose gels to test the amplification success and confirm fragment length. PCR products (single bands) were purified using an enzymatic cleanup protocol (EXOFastAP). 5 µl of the PCR product was mixed with 0.6 µl of alkaline phosphatase FastAP (1U/µl, Thermo Fisher Scientific) and 0.3 µl of Exonuclease I (20 U/µl) on ice, and incubated in a thermocycler (15 min at 37°C and 15 min at 85°C). Purified PCR products were used for bidirectional Sanger sequencing performed by Macrogen Europe. Our laboratory workflow was carefully carried out to minimize potential background contamination from air, dust, laboratory surfaces, equipment and kit reagents. Pre-sterilized disposable materials and equipment were used throughout. All surfaces and non-disposable glass or plexiglass materials were bleach-cleaned (e.g., beakers, pipettes, counting dishes, etc.) and sieves were sonicated for 5 min in between each sample processing. The Tungsten wire was cleaned with 70% ethanol followed by flame sterilization prior to picking each new specimen. To ensure a reliable species identification, the 18S sequences were blasted against Genbank to identify the nematodes using DNA barcoding. Individuals identified as Litoditis marina were further identified to the correct cryptic species using the quantitative real-time PCR (qPCR) method with species-specific primers for the ribosomal internal transcribed spacer region (ITS) (see Derycke et al., 2012). Melting curves were analyzed for the quality control of the results.

2.3 16S rRNA Gene Amplification and Library Preparation for Illumina Miseq Sequencing

We used an adapted protocol of the Earth Microbiome project for the bacterial ribosomal DNA amplification and library preparation for the Illumina Miseq sequencing (Caporaso et al., 2012; Gilbert et al., 2014). DNA extracts of 345 single nematodes were used for the amplification of 16S rRNA gene fragments. The 515F forward (5’-GTGCCAGCMGCCGCGGTAA-3’) and 806R reverse (5’-GGACTACHVGGGTWTCTAAT-3’) PCR primers were used to target the V4 region (ca. 390 bp) (Caporaso et al., 2012). The forward primer contained the 5’ Illumina adapter, the forward primer pad and primer linker, and the 515F primer. The reverse primer contained the reverse complement of the 3’ Illumina adapter, the reverse primer pad and primer linker, the 806R primer and a Golay barcode, unique for each sample (Caporaso et al., 2012), in order to identify each sample during downstream in silico analysis. PCR amplification was performed in 20 µl total reaction volume containing 11.4 µl water, 4 µl buffer (10X), 0.4 µl dNTPs (10mM), 0.2 µl Phusion (high fidelity) polymerase, 1 µl of each primer (10 μM) and 2 µl DNA template. PCR cycle conditions included: an initial denaturation step of 30 s at 98°C, followed by 35 cycles of 10 s at 98°C, 30 s at 65°C, 15 s at 72°C, and a final step of 10 min at 72°C (Derycke et al., 2016). Samples were amplified in technical triplicates, which were then pooled in one sample. Pooled PCR products were purified using a magnetic bead purification protocol with Agencourt AMPure XP beads (Beckman Coulter, CA, USA; the manufacturer’s protocol was followed). PCR library concentrations of the purified samples were measured with a Qubit 3.0 fluorometer (Thermo Fisher Scientific) and equimolarly pooled (15ng). Three libraries-pooled samples were prepared for 3 Illumina Miseq runs, which contained 100, 112 and 133 pooled nematode individuals, respectively. Separate runs were used because of the number of individual samples and the availability of Golay barcodes, necessary for the identification of individual samples during data analysis. Quality control of the three libraries was done with the Bioanalyzer (Agilent Technologies). Purification with AMPure XP beads of the libraries was repeated to remove primer dimers for pools two and three. Each library was sequenced separately on the Illumina Miseq (2 × 300 bp paired-end runs) by Edinburgh Genomics (UK).

2.4 Data Analysis

The demultiplexed Illumina paired-end sequences of the three sequencing runs were further processed together using the DADA2 pipeline (Callahan et al., 2016a) in R, v. 4.0.1 (R Core Team, 2020) and filtered/trimmed using default settings and the Bioconductor workflow (Callahan et al., 2016b), except for truncating of forward and reverse reads at 240 and 200 bp, respectively. Forward and reverse reads with >2 estimated errors were filtered and were truncated at both ends, where read quality dropped below a quality score of 2 (TrunQ = 2). Trimming and filtering was performed on paired reads jointly (both reads must pass the filter for the pair to pass). Adapter sequences were removed, denoised paired-end sequences were merged and Amplicon Sequence Variants (ASVs) were generated after separate estimations of error rates on forward and reverse reads. A Chimera removal step was performed using the consensus method with the ‘removeBimeraDenovo’ function and taxonomy was assigned to ASVs with the Naïve Bayesian RDP classifier (Wang et al., 2007) using the ‘assignTaxonomy’ function with default settings (50% confidence) and the SSU Silva database (v132) as a reference. After filtering and quality control steps, 5 samples were removed based on number of reads and therefore, the full dataset further analyzed included sequences of 340 individual nematode specimens.

2.5 Microbial Diversity Analysis

Downstream data analysis and visualization were conducted using the package ‘phyloseq’ in R (McMurdie and Holmes, 2013). Microbiome diversity was assessed in terms of alpha and beta diversity, which were then statistically compared according to research questions: (i) among species per microhabitat (1 factor: Species); (ii) among microhabitats per species (1 factor: Microhabitat); (iii) across microhabitats and time per species, for species that occurred in the same microhabitats in both sampling seasons (2 factors: Microhabitat, Time); and (iv) among species and microhabitats of different decomposition status (2 factors: Species, Status) for Ulva sp. fresh vs decomposing thalli and for Spartina anglica attached vs detached leaves. For the last research question (iv), two subsets of species co-occurring on these different microhabitats were used and microbiome diversity and composition within microhabitat were compared between ‘Species’ (Sp) and ‘Status’ (St) in a cross-factorial design. Note here that these are in fact different microhabitats per se, but we use the term ‘Status’ here for the purpose of comparison. For the above comparisons, except for (iii), data from the autumn and spring sampling occasions were analyzed separately, due to a lack of sufficient specimens from the same microhabitat on one of the two sampling seasons for most species.

Alpha diversity parameters (i.e., taxon richness as number of observed ASVs, Shannon diversity index) were estimated using the rarefied dataset at a minimum sequencing depth of 10,000 reads (32 samples with lower sequence numbers were excluded). This threshold was chosen as a trade-off to keep as many samples as possible in the dataset while the number of taxa observed in the retained samples was still in the plateau phase of the rarefaction curve (see SM2-Figure 1A). Alpha diversity of the rarefied dataset was then compared between grouping factors using Analysis of Variance (ANOVA), the assumptions of normality of the data and homogeneity of variances, tested using Shapiro–Wilk test and Levene’s test respectively, being met. Microbial community composition analysis was based on the Generalized UniFrac distances (GUniFrac) calculated with the ‘GUniFrac’ package (Chen et al., 2012), using the rarefied dataset. Generalized UniFrac distance was used instead of weighted or unweighted UniFrac distance, containing an extra parameter α (here α = 0.05) to weigh on abundant lineages, so that the distance is not dominated by highly abundant ASVs (Chen et al., 2012). Comparisons of microbiome composition across grouping factors were performed using Permutational Multivariate Analysis of Variance (PERMANOVA) based on the GUniFrac distances, using the ‘vegan’ package (Oksanen et al., 2017). Pairwise difference tests were performed on significant factors, and multivariate dispersions were tested with PERMDISP. Principal coordinates analysis (PCoA) was used to visualize community composition differences within and between grouping factors, using the ‘ade4’ package (Dray and Dufour, 2007). Differential abundance analysis for sequence count data was performed for the complete dataset (non-rarefied sequences of all species and microhabitats) using ‘DESeq2’ with standard parameters (Love et al., 2014) and a significance threshold of a=0.05 (adjusted p-value). Further plotting was performed using the R packages ‘microbiome’ (Lahti and Shetty, 2012-2019) and ‘ggplot2’ (Wickham, 2016) and basic R functions.

3 Results

3.1 Taxonomic Characterization of Microbiomes Associated With Marine Nematodes

The microbiomes of a total of 340 nematode specimens belonging to 10 species (5 genera, 3 families) from 6 different microhabitats (i.e., Fucus spiralis, Fucus vesiculosus, Ulva sp. fresh and decomposing thalli, and Spartina anglica attached and detached leaves) at two sampling moments/seasons (autumn and spring) were characterized (Table 1) and resulted in a total of 28,620,125 filtered reads for the entire dataset with an average sequencing depth of 84,176 reads (median: 68,447 reads) per nematode specimen (see SM2: Table 1). A total of 14,117 distinct ASVs were obtained (average 135 ± 106 per specimen). The vast majority of ASVs were assigned to the Bacterial domain (99.05%, 13,984 ASVs), a few were Archaea (0.34%, 49 ASVs) and the rest remaining taxonomically unassigned.

The phylum Proteobacteria was the most abundant bacterial phylum regardless of the nematode species (46.45% of the bacterial taxa), followed by Bacteroidetes (20.47%), Actinobacteria (5.97%), Firmicutes (5.67%), and Cyanobacteria (2.47%; Figure 1A). Within the phylum Proteobacteria, a high proportion of ASVs was assigned to the Gammaproteobacteria (48.40%), with the families Halomonadaceae (12.15% of Gammaproteobacteria), Shewanellaceae (3.85%) and Pseudomonadaceae (2.58%) and the genera Halomonas, Shewanella and Pseudomonas as the most common representatives (Figure 1B; SM2-Figures 5, 6). The second most abundant class within Proteobacteria was Alphaproteobacteria (39.50%), with the family Rhodobacteraceae being the most abundant within the class (33.01%), followed by Rhizobiaceae (14.81%; Figure 1C; SM2-Figures 5, 6). Within the phylum Actinobacteria, more than 70% of the ASVs belonged to the Actinobacteria class and were mainly assigned to the families Corynebacteriaceae and Micrococcaceae (more than 10 families with low prevalence < 0.05) (Figure 1D). Within the phylum Firmicutes, the most abundant ASVs were assigned to the bacterial class Bacilli (52.33%), more specifically the families Staphylococcaceae (22.41% of Bacilli) and Streptococcaceae (17.35%), represented by the genera Staphylococcus and Streptococcus, respectively (Figure 1E; SM2-Figures 5, 6). The ASVs assigned to the phylum Bacteroidetes belonged to the class Bacteroidia (96.72%), mostly represented by the families Flavobacteriaceae (34.54%), Saprospiraceae (17.11%) and Weeksellaceae (3.77%; Figure 1F; SM2-Figures 5, 6). Archaea were associated with only a few specimens of some nematode species and were mainly represented by the genus Candidatus-Nitrosopumilus, Methanobacterium (only specimens of H. disjuncta Gd3) and Methanobrevibacter (specimens of D. meyli and H. disjuncta Gd3; SM2-Figure 4).

Figure 1 Relative abundance of the most abundant bacterial taxa (prevalence 0.5) of nematode microbiomes per species collected from all six microhabitats and both seasons. ASVs were assigned to phylum for all Bacteria (A) and to family level within the classes Gammaproteobacteria (B) and Alphaproteobacteria (C) and within the phyla Actinobacteria (D), Firmicutes (E) and Bacteroidetes (F). Nematode species names (see Table 1 for full names) are abbreviated as Dd, D. delyi; Dm, D. meyli; Dmd, D. cf. meyli; Eu, Eumonhystera sp.; Gd1/3, H. disjuncta Gd1/3, Papa, P. paetzoldi; Pasp1, Panagrolaimus sp.1; Pm1/2, L. marina Pm1/2. ‘Other’ refers to aggregated taxa with relative abundance < 0.005.

3.2 Microbiome Diversity of Co-Occurring Nematode Species in the Same Microhabitats

3.2.1 Comparing Nematode Microbiomes Among Their Shared Microhabitats

No significant differences were revealed in microbial diversity (Shannon-Wiener diversity) and taxon richness among co-occurring nematode species in each microhabitat, for five microhabitats collected in autumn and four collected in spring, with the only exception of microbial diversity in nematodes obtained from detached decomposing leaves of Spartina in autumn, where Eumonhystera sp. showed significantly lower Shannon diversity than D. meyli (see SM1-Table 1).

Microbiome composition, in contrast, differed significantly among co-occurring species in half of the microhabitats in both seasons (3 out of 6 microhabitats in autumn and 2 out 4 in spring) (Table 2; Figure 2; see SM1-Table 1). However, in most cases, variances were significantly heterogeneous (PERMDISP: p<0.05; Table 2; see SM1-Table 1), hence the observed differences have to be interpreted with caution. Specifically, microbial community composition was significantly different among nematode species co-occurring on decomposing F. spiralis (p = 0.013), F. vesiculosus (p = 0.009) and Ulva sp. (p = 0.003) in autumn (Table 2; Figure 2). In spring, significantly different microbiome composition was revealed for nematode species co-occurring on decomposing F. spiralis (p=0.024) and on fresh Ulva sp. (p=0.006), with significant heterogeneity of dispersions for Fucus but not for Ulva (PERMDISP: p=0.021,p=0.37 resp.; Table 2; Figure 2). The significant pairwise differences observed per microhabitat were not related to specific pairs of species and appeared unrelated to the taxonomic relatedness of the nematode species (for detailed pairwise comparisons, see SM1-Table 2). Differential abundances of microbial genera or families tested only explained 14% or 16%, respectively, of the variation of microbiomes among microhabitats for all the different nematode species together (see SM2-Figure 7).

Table 2 Table of summary statistics of PERMANOVA and PERMDISP comparing microbiome composition of nematode specimens among species in each microhabitat where they co-occurred in autumn and spring.

Figure 2 Principal coordinates analysis (PcoA) plots based on generalized UniFrac distances of the microbiome of single nematode specimens of co-occurring nematode species on Fucus spiralis in (A) autumn and (B) spring, on (C) decomposing Ulva in autumn and (D) fresh Ulva in spring, and on (E) Fucus vesiculosus in autumn. Nematode species are represented by different colours and abbreviated as Dd: D. delyi, Dm: D. meyli, Dmd: D. cf. meyli, Eu: Eumonhystera sp., Gd3: H. disjuncta Gd3, Papa: P. paetzoldi, Pasp1: Panagrolaimus sp.1, Pm1/Pm2: L. marina Pm1/Pm2.

3.2.2 Comparing Nematode Microbiomes Between Microhabitats in Different Decomposition Status

From the above microhabitats, two were additionally examined to compare microbiome differentiation in relation to their status of decomposition. The first subset included two monhysterid nematode species: D. meyli and H. disjuncta Gd1, and a rhabditid species: L. marina Pm2, which were encountered in the field on both fresh and decomposing Ulva sp. in autumn. The second subset included two monhysterid nematode species: D. meyli and Eumonhystera sp., which were encountered on both attached and detached Spartina leaves in autumn. No significant differences in Shannon diversity and ASV richness were revealed across different microhabitat status and species for both test cases (all ANOVAs: p>0.05; SM1-Table 3). Microbiome composition, by contrast, differed significantly between specimens collected from Ulva sp. in different decomposition status (see SM2-Figure 2) regardless of nematode species and between nematode species regardless of the decomposition status (Sp: p=0.024; St: p=0.001; Sp × St: p=0.093; PERMDISP: p>0.05; p>0.05 and p=0.003, respectively; SM1-Table 3; SM2-Figure 2A). Pairwise comparisons were only significant for two species found associated with Ulva: Pm2 vs D. meyli (SM1-Table 3). In contrast, no significant differences in microbiome composition were observed for the Spartina test case (SM1-Table 3).

3.2.3 Comparing Nematode Microbiomes Between Sampling Seasons

When looking at microbiome differences between sampling seasons, no differences were detected in Shannon diversity or microbiome composition for all species tested across microhabitats and/or sampling season (p>0.05 for all), except for one: D. cf. meyli, which exhibited a significant difference in microbiome composition as a function of the interaction between sampling season (S) and microhabitat (M) (M × S: p=0.019; M and S: p>0.05; PERMDISP: p>0.05 for all; SM1-Table 4). Significant pairwise differences were detected between seasons for specimens found on fresh Ulva (pairwise PERMANOVA: p=0.017) but not on Fucus or Spartina (p>0.05; SM1-Table 4). Moreover, significantly higher taxon richness was revealed for H. disjuncta Gd3 in autumn than in spring, regardless of the type of microhabitat (ANOVA: S: p=0.015; M and M × S: p>0.05), but no differences were detected in Shannon diversity or microbiome composition for the species across microhabitats and sampling seasons (S, M and M × S: p>0.05, PERMDISP: p>0.05 for all; SM1-Table 4).

3.2.4 Microbiome Diversity of Nematode Species Between Microhabitats

Out of 10 nematode species, only two showed significant differences in microbiome composition between the microhabitats they inhabited, i.e. L. marina Pm2 among 4 microhabitats in autumn (p=0.001; PERMDISP: p=0.043) and Panagrolaimus sp.1 between detached Spartina leaves vs fresh Ulva thalli in spring (p=0.01; PERMDISP: p=0.651; Figure 3). For L. marina Pm2, significant pairwise differences were revealed between decomposing vs fresh Ulva sp. (p=0.007) and between decomposing Ulva and the two Fucus spp. in autumn (p=0.026 and p=0.002 for F. spiralis and F. vesiculosus, resp.; Figure 3A; see SM1-Table 4). Pm2 also exhibited significant differences in taxon richness between microhabitats (ANOVA:p=0.009), more specifically between decomposing Ulva and F. vesiculosus in autumn (p=0.006; Figures 4A,B) and Panagrolaimus sp.1 exhibited significantly higher microbial diversity on Spartina than on Ulva in spring (ANOVA:p=0.025; Figures 4C,D). No significant differences in microbiome diversity and composition were detected between microhabitats for other species on either sampling season (see SM1-Table 5; SM2-Figure 3). This is likely a consequence of variances within groups being significantly heterogeneous for several species (6 out of 10 in autumn, 1 in spring), indicating a high intraspecific variability of the microbiome composition, regardless of microhabitat. This was the case for D. cf. meyli, Eumonhystera sp., H. disjuncta Gd1 and Gd3, L. marina Pm2 and Panagrolaimus sp.1 in autumn and Eumonhystera sp. in spring (PERMDISP: p<0.05 for all; SM1-Table 5). Differential abundance analysis at the microbial genus level showed significant differences between the aforementioned nematode species pairs, but overall explained only 14% of the variation among individual microbiomes for the complete dataset irrespective of the nematode taxonomic identification (named species here) and the microhabitat type (Figure 5).

Figure 3 Principal coordinates analysis (PcoA) plots based on generalized UniFrac distances of the microbiome of (A) Litoditis marina Pm2 specimens from decomposing thalli of Fucus spiralis (in black), Fucus vesiculosus (in yellow) and from fresh and decomposing Ulva thalli (in green and in blue, resp.) in autumn and of (B) Panagrolaimus sp.1 specimens from detached Spartina leaves (in black) and fresh Ulva thalli (in yellow) in spring.

Figure 4 Box-whisker plots of alpha diversity measures (Shannon diversity index and taxon richness) for (A, B) Litoditis marina Pm2 in autumn and for (C, D) Panagrolaimus sp.1 in spring, grouped per microhabitat where the nematode species occurred (i.e. Fucus spiralis, Fucus vesiculosus, Ulva decomposing, Ulva fresh for Pm2 and detached Spartina leaves and fresh Ulva for Panagrolaimus sp).

Figure 5 Principal coordinates analysis (PCA) plots based on the differential abundance data of the microbiome of the different nematode species (color) and microhabitats (shape), expressed at the taxonomically assigned genus level. Nematode species names (see Table 1 for full names) are abbreviated as Dd, D. delyi; Dm, D. meyli; Dmd, D. cf. meyli; Eu, Eumonhystera sp.; Gd1/3, H. disjuncta Gd1/3; Papa, P. paetzoldi; Pasp1, Panagrolaimus sp.1; Pm1/2, L. marina Pm1/2. Microhabitats are abbreviated as Fucus_spira, Fucus_vesi for Fucus spiralis and Fucus vesiculosus, Spart_deco, Spart_old for detached and attached Spartina leaves and Ulva_deco, Ulva_fresh for decomposing and fresh Ulva, respectively.

4 Discussion

4.1 The Microbiome of Bacterivorous Nematodes Is Highly Diverse

Our data revealed a highly diverse microbial community (13,984 ASVs were assigned to bacterial taxa) associated with intertidal nematode specimens of different taxonomical identities (average min.120 and average max. 197 ASVs were assigned per specimen across 10 nematode species) and a high number of reads (see SM2-Figure 1). The average number of raw reads per sample (84,176 reads) was high, although comparisons with other studies on microbiome of nematode individuals from marine habitats can not be directly made due to different methods (e.g. average 31,790 reads using the Illumina Miseq platform and OTUs in Schuelke et al. (2018); 5,744 reads using the 454 GS-FLX platform and OTUs in Derycke et al., 2016). Microbiomes of bacterivorous nematodes comprise not only bacteria that are part of their diet (microbiome ‘sensu lato’), but also host-associated bacteria that hold ecological interactions (e.g., commensal, symbiotic or pathogenic) with nematodes (microbiome ‘sensu stricto’) (Musat et al., 2007; Cabreiro and Gems, 2013; Derycke et al., 2016). Gut-symbiotic bacteria (i.e., gut microflora) can also be related to resource utilization of bacterivorous nematodes and can substantially contribute to their microbiomes (Baquiran et al., 2013; Cabreiro and Gems, 2013). Hence, we assume that a significant part of the microbiome of our nematode specimens are either food-derived or potentially involved in the resource utilization process, as also suggested by Derycke et al. (2016).

4.2 Microbiomes of Co-Occurring Bacterivorous Nematode Species Display Strong Overlap and a High Intraspecific Variability

Niche differentiation is a mechanism enabling local coexistence of closely related species (Chesson, 2000; Leibold and McPeek, 2006). We anticipated a degree of microbiome differentiation among species that share the same microhabitat, resulting at least in part from a differential resource use. Our results do not provide strong support for this hypothesis, as microbiomes did not show a significantly distinct composition among co-occurring species in half the number of microhabitats examined in autumn and in spring, exceptions being the two Fucus species and decomposing Ulva sp. in autumn, and F. spiralis and fresh Ulva sp. in spring.

Nevertheless, differences in microbiome composition of nematode species related with different stages of habitat decomposition were observed in our study. The state of decomposition per se means that some bacterial activity is already in place and thus slight differences in microbiomes for some nematode species could be expected. Yet, differences in nematode microbiomes related to decomposition status were only observed for species found on green macroalgae Ulva sp. Different microhabitat preferences related to different stages in the decomposition process have been suggested for closely related bacterivorous species associated with cordgrass or macroalgae. Specifically, Diplolaimelloides meyli was suggested to be associated with early-stage, and Diplolaimella dievengatensis and Diplolaimelloides oschei with later stages of the decomposition process of Spartina leaves (Moens et al., 1999). Litoditis marina is known as a particularly successful colonizer of early-stage decomposition of macroalgae (Derycke et al., 2007b) and Pm2 exhibited differences in microbiome composition between early and later decomposition stage of Ulva in our study. That was also the case for D. meyli and Halomonhystera disjuncta - Gd1, suggesting that resource divergence may be related to some extent to microhabitat partitioning.

In agreement with previous studies on L. marina and C. elegans (Derycke et al., 2016; Dirksen et al., 2016), we observed a high intraspecific variability in microbiome composition, which may provide an explanation for the lack of clear patterns of microbiome divergence among species. When interpreted as a consequence of differences in the microbiome ‘sensu lato’, this high intraspecific variability in microbiome composition may result from a substantial intraspecific competition (Svanbäck and Bolnick, 2007). Strong intraspecific competition may prevent individuals from exploiting the whole range of available resources (Svanbäck and Bolnick, 2007; Guden et al., 2021). The resulting interindividual resource specialization can expand the niche breadth of the total population (Bolnick et al., 2007; Svanbäck and Bolnick, 2007). Population niche breadth is thought to represent a balance between the diversifying effect of intraspecific competition and the constraints imposed by interspecific competition. Trade-offs between resource utilization and competitive success have been suggested to favor generalist populations with specialist individuals (Bolnick et al., 2003). From that perspective, the combination of a large interindividual microbiome variation with limited interspecific microbiome divergence may be suggestive of a general importance of intraspecific competition in the communities of bacterivorous nematodes studied here.

Species-specific nematode-microbiota associations may allow organisms to expand their resource utilization or niche. It has been suggested that putative mutualists may facilitate food digestion by nematodes through secretion of toxins, antibiotics or digestive enzymes, but they may also protect their hosts against pathogenic bacteria (Baquiran et al., 2013). Microbiome taxa may also be involved in mediating differential tolerances of their hosts to environmental conditions (Derycke et al., 2016; Dirksen et al., 2016), which can also be important for niche partitioning. For instance, bacteria of the genus Acinetobacter were found in association with Pm1 (and Pm3) by Derycke et al. (2016) and in the present study, have been suspected to be involved in abiotic tolerances of Pm1 (Derycke et al., 2016). Putative symbiotic associations of free-living marine nematodes and specific bacterial/archaeal taxa have also been found related to important ecological functions (e.g. oxidizing, nutrient cycling), which suggests that although these relationships remain largely unexplored, host-associated microbiomes can diversify species niche (Schuelke et al., 2018). Therefore, it is possible that symbionts play a role in niche partitioning of co-occurring nematodes, not only through their contribution to their host’s resource use, but also through influences on other aspects of their niche.

4.3 The Microbiome of Most Nematode Species Does Not Differ Between Microhabitats

The microbiomes of most nematode species did not differ among microhabitats, in contrast with our expectations of microbiome divergence as an indication of differential resource availability. On the contrary, each species showed a large overlap of microbiome diversity and composition across the microhabitats they occupied. There were only two (Rhabditida) nematode species in our study with clear differences in their microbiome depending on the microhabitat from which they were collected. The lack of nematode microbiome differentiation among microhabitats could be explained by two hypotheses: (a) a lack of microbiome differences associated with the microhabitats, and (b) the high intraspecific variability in microbiome composition of the nematodes. The first hypothesis is less probable as clear differences in epibiotic microbial communities associated with species and even (sub)structures of macroalgae have been reported in the literature (Lachnit et al., 2009; Lachnit et al., 2011), demonstrating algal species-specific and structure-specific microbiomes. Those microbiome differences between microhabitats, however, are not reflected in the microbiomes of the nematodes, which suggests that the nematode species studied here either utilize the same resources in different microhabitats, or that their resource use is so variable that no significant differences occur. Another study on the microbiomes of a large number of free-living marine nematode species showed a lack of habitat-microbiome correlation not only for bacterivorous nematodes but also for other functional groups (i.e. herbivores, predators, omnivores; Schuelke et al., 2018). According to their results, the structure and diversity of microbial communities associated with free-living marine nematodes from various geographic locations were not correlated with habitat type, albeit they were collected at much larger spatial scales (e.g., intertidal sediments, methane seeps, deep sea). Our data extent these results and show that also nematode species belonging to the same functional group (bacterivores) and collected from different habitats in one geographic location have microbiomes that are unlinked to (micro)habitat type. We therefore agree with the suggestion by Schuelke et al. (2018) that free-living nematodes in marine habitats may rather exhibit microbiome composition plasticity and be generalists in their feeding, which may also be the case for the bacterivorous nematodes examined in our study. Studies on entomophile nematodes have also demonstrated a lack of correlation of nematode-associated microbiomes and their respective (micro)habitats (Koneru et al., 2016; Meyer et al., 2017). By contrast, diverse microbiomes of plant-parasitic nematodes and of the well-studied soil nematode C. elegans did differ in relation to their microbial environments (Berg et al., 2016; Elhady et al., 2017). Microbiome structure was also related to habitat (and species) for freshwater and soil tardigrades and marine isopods (Mattila et al., 2014; Vecchi et al., 2018). On the other hand, the high intraspecific variation observed in microbiome composition suggests that nematodes more likely do not utilize the same resources but extend their trophic niche by individual resource specialization, which explains the lack of significant microbiome differentiation between microhabitats. Future research on nematode microbiomes, including a broader range of habitats and other feeding guilds (i.e. detritivores, predators, omnivores), would be interesting to draw clearer conclusions on the role of trophic niche specialization of nematodes for their coexistence patterns.

Although we could test temporal variability in nematode microbiomes for only four species, all belonging to the same family (Monhysteridae), there was only one observed case of temporal variability in resource use which may result in temporal niche differentiation among co-occurring species on Ulva sp. Priority effects may also play a role during colonization of macroalgae (Derycke et al., 2007b) or any substratum by competing species, which can also result in temporal differentiation in resource utilization. A pattern of temporal differentiation in microhabitat colonization has also been reported for co-occurring Litoditis marina cryptic species in the field (Guden et al., 2018).

Conclusion

Co-occurring free-living nematodes from an intertidal habitat did not show clear resource differentiation patterns among species sharing the same microhabitats, nor did single species among the different microhabitats they occupied. Contrary to our expectations, microbiomes of bacterial-feeding nematodes were not strongly species-specific nor microhabitat-specific, apart from an indication that nematode microbiomes differ in relation to the decomposition process of their phytal substratum. However, a pronounced intraspecific variability in microbiome diversity and composition was observed for the examined estuarine nematode species. We suggest that this intraspecific variability, most likely resulting from intraspecific competition, is potentially the main driver promoting coexistence of sympatric species.

Data Availability Statement

The data that support the findings of this study are publicly available online: Sanger 18S sequences at BOLD dataset [DS-MOLECOL] and GenBank [accession numbers: ON713487 - ON713781] and HTS raw reads at NIH Sequence Read Archive [BioProject “Microbiome partitioning of coexisting nematode species” - accession number: PRJNA835312].

Author Contributions

A-MV designed and performed this research, the data analyses and wrote the paper. TM and SD contributed to the design of the research and the writing of the manuscript. AR contributed to lab analyses, NM to data analyses and RG together with all co-authors revised the manuscript. TM provided funding for this research. All authors contributed to the article and approved the submitted version.

Funding

Financial support for this study was obtained from the Flemish Science Fund FWO through projects G0C3218N and 3G.0192.09, and from the Special Research Fund of Ghent University (BOF) through project 17GOA-026. The research leading to results presented in this publication was carried out with infrastructure funded by FWO (project G0H3817N) in the framework of EMBRC Belgium - international research infrastructure I001621N.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

High-throughput sequencing was carried out by Edinburgh Genomics, The University of Edinburgh. This article is based on work that was part of the Ph.D. thesis of AMV.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2022.881566/full#supplementary-material

References