The marine environmental microbiome mediates physiological outcomes in host nematodes

Abstract

Background

Nematodes are the most abundant metazoans in marine sediments, many of which are bacterivores; however, how habitat bacteria affect physiological outcomes in marine nematodes remains largely unknown.

Results

Here, we used a Litoditis marina inbred line to assess how native bacteria modulate host nematode physiology. We characterized seasonal dynamic bacterial compositions in L. marina habitats and examined the impacts of 448 habitat bacteria isolates on L. marina development, then focused on HQbiome with 73 native bacteria, of which we generated 72 whole genomes sequences. Unexpectedly, we found that the effects of marine native bacteria on the development of L. marina and its terrestrial relative Caenorhabditis elegans were significantly positively correlated. Next, we reconstructed bacterial metabolic networks and identified several bacterial metabolic pathways positively correlated with L. marina development (e.g., ubiquinol and heme b biosynthesis), while pyridoxal 5’-phosphate biosynthesis pathway was negatively associated. Through single metabolite supplementation, we verified CoQ₁₀, heme b, acetyl-CoA, and acetaldehyde promoted L. marina development, while vitamin B6 attenuated growth. Notably, we found that only four development correlated metabolic pathways were shared between L. marina and C. elegans. Furthermore, we identified two bacterial metabolic pathways correlated with L. marina lifespan, while a distinct one in C. elegans. Strikingly, we found that glycerol supplementation significantly extended L. marina but not C. elegans longevity. Moreover, we comparatively demonstrated the distinct gut microbiota characteristics and their effects on L. marina and C. elegans physiology.

Conclusions

Given that both bacteria and marine nematodes are dominant taxa in sedimentary ecosystems, the resource presented here will provide novel insights to identify mechanisms underpinning how habitat bacteria affect nematode biology in a more natural context. Our integrative approach will provide a microbe–nematodes framework for microbiome mediated effects on host animal fitness.

Background

Multicellular organisms, including animals and plants, interact closely with diverse and abundant microbial communities in their natural habitats [1,2,3]. Both the microbes physically associated with the host and the animal gut microbiota influence host physiology and fitness [2,3,4]. The microbiota communities modulate essential biological processes such as growth, development, reproduction, longevity, and responses to biotic and abiotic stresses, in their host animals and plants [5,6,7]. Changes in environmental bacteria can affect the behaviors of zebrafish such as locomotion [8] and anxiety-related behavior [9]. Experimental studies in mouse models have demonstrated substantial impacts of the gut microbiota on neuron development, angiogenesis, immunity, and disease [10, 11]. The dysbiosis of gut microbiota has been associated with a wide spectrum of human diseases such as obesity, diabetes, cancer, and neurological disorders [12].

Among the model organisms utilized to study the interactions between microbiota and host animals, the terrestrial nematode Caenorhabditis elegans is an excellent model to study microbiome-mediated functions, due to its unique characteristics such as a short generation time, clear genetic background, and advanced functional genomics resources [13]. C. elegans encounters different bacteria in its native habitats, which could be potentially used as dietary resources, some of which colonize the worm’s gut [14]. It was reported that nearly 80% of the > 550 habitat bacterial isolates could individually support C. elegans development [15]. The effects of several microbial metabolites on host physiology have been investigated in C. elegans, e.g., vitamins [16,17,18], siderophores [19, 20], and neurotransmitters [21]. Genome-scale metabolic models were recently reconstructed for 77 members of C. elegans natural habitat bacteria, revealing key bacterial metabolic traits related to gut bacterial colonization and host fitness [22]. Recently, the gut microbiota of C. elegans was reported to be shaped by both host genetics and environmental factors [14, 22,23,24].

Although both nematodes and bacteria are dominant organisms in marine sedimentary ecosystems, and many marine nematodes are bacterivores, however how native bacteria affect the essential biological processes of marine nematodes is largely unexplored. Litoditis marina, a marine free-living bacterivore nematode, is widely distributed and plays an important role in marine ecosystems [25,26,27]. We have recently established L. marina as a promising marine model organism, including a short generation time, easy to be cultured in the lab, multiple inbred lines with clear genetic background, high-quality genome assembly, annotations, and functional genomics resources [26, 27]. It was reported that the gut microbiome of cryptic species of L. marina is highly diverse and species-specific [25]. Native bacteria might provide nutrient resources to L. marina, and the natural microbiome is affected by global climate change, however, the composition of the L. marina habitat microbial communities and the natural microbiome-mediated functions on nematode fitness are largely unknown.

In this study, we characterized the native bacterial community compositions and isolated 539 bacterial strains from L. marina natural habitats. Next, the impacts of 448 single native bacterial isolates on nematode development were assessed, and a representative set (73 bacteria, termed HQbiome) was further examined for their effects on host physiology, of which 72 whole genome sequences were generated. Unexpectedly, we found that the effects of marine native bacteria on the development of L. marina and its terrestrial relative C. elegans were significantly positively correlated. Based on the whole genome sequences of HQbiome, we reconstructed bacterial metabolic networks and assessed their metabolic capacities. Based on HQbiome bacterial metabolic networks reconstruction with random forest regression analysis, together with single metabolite supplementation, we demonstrated that CoQ10, heme b, acetyl-CoA, and acetaldehyde promoted L. marina development, while vitamin B6 attenuated growth. Notably, we found that only four growth development correlated metabolic pathways were shared between L. marina and C. elegans. Strikingly, we found that glycerol supplementation significantly extended L. marina but not C. elegans longevity. Moreover, we demonstrated the gut microbiota characteristics and bacteria affect the physiology of both L. marina and C. elegans.

Results

The natural bacterial microenvironment of L. marina

To characterize the bacterial community composition that L. marina encounters in its natural environment, we performed full-length 16S rRNA gene amplicon sequencing of 80 intertidal sediment samples from the L. marina habitats (Huiquan Bay, Qingdao, China) over a 13-month period (Additional file 1: Fig. S1). A total of 5783 operational taxonomic units (OTUs) were identified, belonging to 3571 unique bacterial species (Additional file 2: Table S1A). Of these, the dominant phyla were Proteobacteria (46.33%), Bacteroidetes (22.03%), Cyanobacteria (11.11%), Actinobacteria (6.30%), Verrucomicrobia (3.13%), Acidobacteria (2.96%), Planctomycetes (1.78%), and Gemmatimonadetes (1.60%) (Fig. 1a and Additional file 1: Fig. S2; Additional file 2: Table S1B, C). Notably, Proteobacteria and Bacteroidetes are also the most abundant phylum under C. elegans natural habitats [15]. As the most abundant phylum, Proteobacteria were mostly enriched with classes such as Gammaproteobacteria (28.18%), Deltaproteobacteria (10.91%), Alphaproteobacteria (5.93%), and Betaproteobacteria (0.87%). We found that phylum Cyanobacteria was significantly abundant in winter (Fig. 1a), which was previously reported to be well adapted to cold environments [28, 29]. Of the 1821 identified genera, the most abundant genus was Woeseia, presenting in all 80 samples, and being the only genus with more than 10% mean abundance (11.80%; Fig. 1b; Additional file 2: Table S1D). Moreover, assessments of α-diversity revealed a significantly lower richness and diversity of the natural microbiome in the coldest month (January) compared with that of other seasons (Additional file 1: Fig. S3; PFDR < 0.05, Kruskal‒Wallis test, Wilcoxon test). Meanwhile, there was large seasonality in the less abundant bacteria (Fig. 1a and Additional file 1: Fig. S3).

[IMAGE OMITTED: SEE PDF]

Bacterial culture collections from L. marina habitats

To ask how individual habitat bacteria may affect the development and longevity of L. marina, we isolated bacteria on simple marine bacterial growth media from sediment samples where L. marina was collected. A total of 3822 colony-forming units (CFUs) were obtained and taxonomically characterized by sequencing the bacterial 16S rRNA gene (Additional file 3: Table S2A), resulting in a total of 539 unique CFUs (Additional file 3: Table S2B), belonging to Gammaproteobacteria (67.06%), Alphaproteobacteria (7.74%), Betaproteobacteria (1.28%), Bacteroidetes (13.34%), Firmicutes (8.84%), and Actinobacteria (1.49%) (Fig. 1c, d). These bacteria isolates were categorized into 10 classes, 34 orders, 68 families, and 192 genera (Additional file 3: Table S2B), among which, all phyla, seven classes, 24 orders, 45 families, 85 genera, and 51 species overlapped with our cultivation-independent full-length 16S rRNA community profiling, respectively (Additional file 3: Table S2C). The majority of 539 CFUs were members of the genera Vibrio, Photobacterium, Pseudoalteromonas, Shewanella, Olleya, Bacillus, and Polaribacter, most of which were also reported from other marine bacterial cultures, which might represent the dominant cultivable bacterial genera in marine and coastal ecosystems [30]. Notably, our cultivation-based collection expanded the diversity of the existing collection, and 28 of them were candidates for new species which had a 16S rRNA gene sequence with ≤ 97% identity to their closest references (Additional file 3: Table S2B).

Impact of environmental bacteria on the development of L. marina

To ask how environmental bacteria affect marine nematode development, we fed L. marina individual marine bacterial isolates (448 of 539 isolates examined in this project). By applying k-means clustering, 121 bacterial strains were categorized as generally “beneficial” (fast growth, percentage of stage 4 larvae (L4) in day 5: ≥ 45.71%), 113 were categorized as “intermediate” (moderate growth, percentage of L4 in day 5: 18.57% ~ 45.71%), and 214 were “detrimental” (slow growth or active killing, percentage of L4 in day 5: < 18.57%) (Fig. 2, Additional file 1: Fig. S4, and Additional file 1: Fig. S5; Additional file 3: Table S2D). Nearly twice as many strains were classified as “detrimental” (47.66%) compared to “beneficial” (27.17%) (Fig. 2a). Specific taxa tend to have predominantly beneficial or detrimental impacts on L. marina, for instance, most Alphaproteobacteria (e.g., Rhodobacteraceae) and Betaproteobacteria (e.g., Comamonadaceae) strains were beneficial to L. marina, while Actinobacteria (e.g., Micrococcaceae and Nocardiaceae), Firmicutes (e.g., Planococcaceae and Bacillaceae), Bacteroidetes (e.g., Flavobacteriaceae), and several Gammaproteobacteria (e.g., Moraxellaceae, Psychromonadaceae, and Vibrionaceae) showed detrimental impacts (Additional file 1: Fig. S6). Note that many of these “detrimental” taxa such as Microbacteriaceae and Colwelliaceae are pathogens of C. elegans and other animals [15]. We found a significant phylogenetic signal for the effects on nematode growth (Pagel’s λ = 0.411, P = 3 × 10−13, Abouheif’s Cmean = 0.192, P = 0.001. The observed phylogenetic signal was robust also to subsampling. Additional file 1: Fig. S7; Additional file 3: Table S2E). However, the distribution of the impacts on L. marina development within the certain phylogenetic tree of our culture collection revealed discrepancies between closely related strains, i.e., different strains belonging to the same family exhibited variant effects, for example, members in Vibrionaceae, Shewanellaceae, Rhodobacteraceae, and Flavobacteriaceae (Additional file 1: Fig. S6).

[IMAGE OMITTED: SEE PDF]

Impacts of marine bacteria on C. elegans egg-laying time

To ask whether marine bacteria effects similar physiological outcomes between marine nematode L. marina and its terrestrial relative, the well-known biomedical model nematode C. elegans, we examined the effects of 448 marine habitat isolates on C. elegans egg-laying time, as a proxy for animal developmental rate. By applying k-means clustering, among 448 bacterial strains, 321 were categorized as “beneficial” (fast growth, egg-laying time: ≤ 73.67 h), 50 were categorized as “intermediate” (moderate growth, egg-laying time: 73.67 ~ 123 h), and 78 were “detrimental” (slow growth or active killing, egg-laying time: > 123 h) (Additional file 1: Fig. S8; Additional file 3: Table S2D). Interestingly, C. elegans likely could utilize these bacteria better compared with marine nematode L. marina, as C. elegans exhibited a higher percentage of beneficial strains (71.65% in C. elegans compared to 27.01% in L. marina). Category “beneficial” include the family Moraxellaceae, Bacillaceae, and Pseudomonadaceae, which is similar (i.e., Moraxellaceae) or different (i.e., Bacillaceae and Pseudomonadaceae) compared to that of L. marina (Additional file 1: Fig. S8; Additional file 3: Table S2D). Similar to L. marina, a weak phylogenetic signal was also observed for the effects on C. elegans egg-laying time (Pagel’s λ = 0.694, P = 1 × 10−15, Abouheif’s Cmean = 0.304, P = 0.001; Additional file 1: Fig. S7c, d; Additional file 3: Table S2E).

Comparing both nematodes revealed that development rates were positively correlated when all 448 habitat isolates are considered (Spearman’s ρ = 0.39, P = 2.2 × 10−16; Additional file 1: Fig. S9). The positive correlation means that strains with beneficial/detrimental impacts on L. marina may also exhibit similar effects on C. elegans, e.g., Gammaproteobacteria, and Actinobacteria, while most Firmicutes strains, which were detrimental for L. marina, also delayed C. elegans egg-laying time. The detrimental impacts of Agarivorans, Colwelliaceae, and Psychromonadaceae in Gammaproteobacteria were particularly pronounced (Additional file 1: Fig. S8b). Notably, the correlation between the impacts on L. marina and C. elegans was not evenly distributed across taxa (Additional file 1: Fig. S9c). For example, some strains belonging to Bacteroidetes, Alphaproteobacteria, and Betaproteobacteria that were beneficial to C. elegans also promoted the developmental rate of L. marina, e.g., Maribacter orientalis, Sulfitobacter faviae, and Comamonas aquatica. However, other members of the same bacterial taxa were detrimental to L. marina (e.g., Lacinutrix undariae, Pacificibacter aestuarii, and Hydrogenophaga taeniospiralis).

Establishment of a representative habitat microbiome HQbiome

We next focused on 72 taxonomically and functionally representative of the abovementioned 448 isolates, together with W. oceani, which was not isolated by our cultivation-based screening, while it is highly represented by our cultivation-independent survey (Fig. 1b), these 73 bacteria were termed HQbiome (Fig. 2b, c; Additional file 3: Table S2F and Additional file 1: Table S5). The representative habitat microbiome HQbiome, covered a range of impacts on L. marina, from beneficial to detrimental, and the proportions of the three categories in HQbiome were similar to that of 448 tested isolates in this study (Fig. 2a, b and Additional file 1: Fig. S5). The integration of the prediction of the metabolic potential of habitat communities by PICRUSt2 and the metabolic model of HQbiome by gapseq revealed that both metabolic competences were positively correlated (Methods; Spearman’s ρ = 0.6, P < 2.2 × 10−16; Additional file 1: Fig. S10a), suggesting that HQbiome could be functionally representative of L. marina habitat bacterial communities. Furthermore, HQbiome strains assembled into communities resembling the habitat microbiome via the cultivation-independent approach and the total cultivation-dependent collection at the phylum level (Fig. 1d and Additional file 1: Fig. S10b).

Of the 73 strains in HQbiome, 29 (39.73%) were “beneficial” to L. marina, including Alphaproteobacteria and Gammaproteobacteria, while 25 strains (34.25%) belong to the “detrimental” category, most of which were Actinobacteria, Bacteroidetes, and Firmicutes (Fig. 2b). The remaining 19 strains were “intermediate”, belonging to Gammaproteobacteria and Firmicutes. Notably, over 61% of L. marina fed Alphaproteobacteria reached L4 on day 5, whereas Actinobacteria and Bacteroidetes strains significantly attenuated animal growth (Fig. 2b). As for the impacts on C. elegans, among 73 HQbiome bacterial isolates, 50 were categorized as “beneficial”, 6 as “intermediate”, and 17 as “detrimental”. Comparing both nematodes revealed that development rates were positively correlated when HQbiome isolates are considered (Spearman’s ρ = 0.57, P = 1.2 × 10−7; Additional file 1: Fig. S11).

Principal coordinates analysis (PCoA) of distances of the impact on the development of L. marina and C. elegans revealed a clear clustering of isolates on the basis of their Gram staining trends, which separated along the first coordinate axis (Additional file 1: Fig. S12c). This separation was mainly explained by the L4 proportion of L. marina on day 5 (Spearman’s ρ = 0.9546475, P = 2.2 × 10−16), which was used as a measure of the developmental rate of L. marina in the present study. Most Gram-positive strains, especially Actinobacteria isolates, impaired L. marina and C. elegans physiology (Additional file 1: Fig. S12a, b). We also investigated the diversity of impacts of HQbiome strains on nematode’s growth within each bacterial phylum (class level for Proteobacteria) to identify bacterial taxa with varying degrees of diverse bacterial groups on nematode physiology (Additional file 1: Fig. S12d). Alphaproteobacteria showed significantly lower diversity compared to Actinobacteria, Gammaproteobacteria, and especially Firmicutes, which exhibit a higher degree of within-taxa diversification of impacts on both L. marina and C. elegans, even though all groups have a comparable degree of phylogenetic relatedness.

Whole genome sequencing and genome-scale metabolic networks reconstruction and validation

We generated whole genome sequences (WGS) for 72 bacterial isolates from HQbiome. Of which, 15 were sequenced using both Oxford Nanopore Technologies (ONT) and Illumina technologies, one with both Pacific Biosciences (PacBio) and Illumina technologies, resulting in either a single-circular chromosome (9 strains) or a single chromosome with 1 to 6 plasmids (Additional file 3: Table S2F). The remaining 56 strains were sequenced with Illumina only, resulting in assemblies with 8 to 211 contigs per genome.

For five genera and two families (Vibrio, Shewanella, Pseudoalteromonas, Pseudomonas, Psychrobacter, Bacillaceae, and Rhodobacteraceae), each including more than three isolates, we identified substantial inter-specific genome variation (Additional file 1: Fig. S13). Based on HQbiome WGS, we constructed a phylogenetic tree and identified unidentical consensus phylogeny compared to the phylogenetic tree inferred from the full-length 16S rRNA gene (Fig. 3a and Additional file 1: Fig. S14). Especially, isolates belonging to Bacillaceae clustered well in the WGS phylogenetic tree.

[IMAGE OMITTED: SEE PDF]

To investigate the gene contents and functions of HQbiome isolates, we annotated their genomes using the COG database [31]. The overall eggNOG-derived functional profile of all HQbiome strains indicated that over 20% of annotated proteins have unknown functions (Additional file 1: Fig. S15a), stressing the need to perform in-depth functional studies. To deeper explore metabolic capacities within HQbiome, we reconstructed genome-scale metabolic models (GEMs) using the gapseq pipeline (Additional file 4: Table S3A). To validate the GEMs, we experimentally quantified the ability of Escherichia coli OP50 and four Vibrio strains to utilize 71 carbon sources using the BIOLOG assay, resulting in 73.20% overlap with the curated models (Additional file 1: Fig. S16; Additional file 4: Table S3M), which was similar to the previous study [22]. In addition, we evaluated the quality of the GEMs by testing if the models can recapitulate carbon source utilization as indicated in the ProTraits database and found that the correct prediction rate is about 90% (37/41), highlighting the high quality of our reconstructed metabolic models (Additional file 1: Fig. S17; Additional file 4: Table S3M). Moreover, quality assessment tests with all GEMs were performed with MEMOTE [32], yielding overall high consistencies: on average, 100% stoichiometric consistency, 99.92% reactions mass-balanced, 99.97% reactions charge-balanced, and 100% metabolite connectivity.

The HQbiome metabolic reconstructions harbored an average of 730 ± 168 genes, 1935 ± 210 reactions, and 1710 ± 147 metabolites (Additional file 1: Fig. S15b). We found that Gammaproteobacteria had larger genome sizes than the other taxa, and Firmicutes had larger sets of genes, reactions, and metabolites than other taxa. Meanwhile, simulation of in silico growth suggests that L. marina developmental rates fed Firmicutes strains were significantly faster than the isolates belonging to other taxa except for Gammaproteobacteria (Additional file 1: Fig. S15b), implying that higher developmental rate was associated with more metabolic reactions and metabolites.

Metabolic capacity and diversity of HQbiome strains

Hierarchical clustering of the genome-informed metabolic potential reflected extensive variability across the HQbiome isolates (Fig. 3b; Additional file 4: Table S3C), exhibiting three main clusters: one with generally “detrimental” effects on L. marina development, comprising Actinobacteria, Bacteroidetes, and Firmicutes (left), and two with generally “beneficial” effects on L. marina growth, including Alphaproteobacteria and Gammaproteobacteria (middle and right). In general, we observed that the basic metabolites essential for their host animals, including nucleotides, essential amino acids, and most cofactors, were found among the “beneficial” isolates (Fig. 3b). Strikingly, the bacterial cluster with detrimental effects on L. marina development lacked the metabolic potential of several lipids (e.g., oleate, stearate, palmitoleate, cyclopropane fatty acid, and 5Z-dodec-5-enoate). The attenuated L. marina development induced by isolates from the “detrimental” cluster could be due to the lack of the abovementioned metabolic capability for lipids and cofactors synthesis.

Moreover, the inferred metabolic competences are consistent with the lifestyle of the L. marina host in the intertidal zone; for example, most HQbiome strains had enzymes that enable tolerance to microaerobic conditions (e.g., cytochrome bd oxidase). Only 50 pathways were common to all 73 HQbiome GEMs (Additional file 4: Table S3A). Metabolic pathways enrichment was detected at different taxonomic levels. For instance, only Shewanella showed the capacity to produce polyunsaturated fatty acid (PUFA; eicosapentaenoate, and docosahexaenoate) and olefins (Additional file 4: Table S3B). Of note, the pathways related to anaerobic respiration (i.e., PWY0-1576, PWY0-1577, and PWY-4601) were particularly enriched in the genera Shewanella. Bacillaceae isolates were able to produce terpenoid (e.g., dammara-20,24-diene and diploterol) and folate (e.g., 6-hydroxymethyl-dihydropterin diphosphate) (Additional file 4: Table S3B). Rhodobacteraceae strains showed metabolic competences for phospholipid synthesis (e.g., phosphatidylcholine, phosphatidylserine, and phosphatidylethanolamine) (three isolates belonging to Pseudomonas and Marinomonas also possess these competences), but apparently Rhodobacteraceae strains lack phenylalanine, serine, and threonine biosynthetic pathways (Fig. 3b; Additional file 4: Table S3B). Moreover, we found that the metabolites with significant presence and absence were observed in certain taxonomic groups and the clusters of effects on nematode growth (Fig. 3b). Actinobacteria strains (e.g., Galactobacter valiniphilus and Arenivirga flava) tended to have reduced genomes (Additional file 1: Fig. S15b) and consequently reduced metabolic capabilities (e.g., lipids), in agreement with the previous report that genome size was an important predictor of the overall biosynthetic capabilities of an organism [33]. Strains belonging to Vibrionaceae and Shewanellaceae lacked arginine biosynthetic pathway, while Rhodobacteraceae lacked phenylalanine, serine, and threonine biosynthetic pathway (Fig. 3b). As expected, lipopolysaccharide biosynthesis capacity was found only in Gram-negative bacteria, whereas teichoic acids and peptidoglycan were found only in Gram-positive strains, which were mainly distributed in the “detrimental” cluster (Fig. 3b, left). Additionally, we identified the presence of specific substrate transporters such as amino acids and cofactors in the HQbiome, and found that several Alphaproteobacteria strains lacked lactate transporters (Additional file 1: Fig. S18).

To further explore the potential behavior of HQbiome isolates in an ecological context, we assessed their genomic and metabolic traits in light of the universal adaptive strategy theory [22] and demonstrated that the stress-tolerating strategy was the most prevalent (Additional file 1: Fig. S15c). There were no significant differences between the three groups (beneficial, detrimental, and intermediate) (Additional file 1: Fig. S19a). Meanwhile, none of the strategies between bacterial supporting nematode development was significant (P = 0.552 and 0.7, Kruskal–Wallis; Additional file 1: Fig. S19b). We subsequently examined the correlation between features of GEMs and host growth phenotypes using Spearman rank correlation (Additional file 1: Fig. S15d). Interestingly, isolates with more exopolysaccharides were significantly positively correlated with L. marina developmental rate, while siderophores were likely detrimental to L. marina.

PCoA of metabolic distances revealed a clear clustering of genomes on the basis of their taxonomy, with distinct clusters for the Proteobacteria and other phyla (Additional file 1: Fig. S20a, PCoA1). Gram-positive and Gram-negative strains were clearly separated, consistent with the phenotypic groupings described above (Additional file 1: Fig. S12c). We then examined the metabolic diversity within each phylum (Additional file 1: Fig. S20b) to identify bacterial taxa with varying degrees of functional versatility. Isolates from Alphaproteobacteria and Actinobacteria showed a significantly lower metabolic diversity compared to those belonging to Firmicutes and Gammaproteobacteria (Additional file 1: Fig. S20b). We subsequently tested the relationship between metabolic and phylogenetic similarities and the bacteria’s metabolic potential, and found that taxonomically related bacteria shared more reactions than taxonomically distant bacteria (Additional file 1: Fig. S20c; Spearman rank correlation, ρ = 0.67, P < 2.2 × 10−16). Phylogenetic similarities appeared to be larger than metabolic similarities, suggesting variations in metabolic competences within taxa (Additional file 1: Fig. S20d). For example, HQB-179 was clearly separated from the cluster of other Alphaproteobacteria isolates, and similar separation was also observed at genus and family levels (HQB-226 clearly separated from HQB-476).

To examine whether the known functional diversity of the HQbiome isolates was captured in our models, we performed in silico growth simulations for the uptake of 55 different carbon sources and the secretion of 41 fermentation products using flux variability analysis (FVA) (Additional file 1: Fig. S21). Almost all isolates could utilize and likely would compete for simple sugars, such as glucose and fructose, while only a few (e.g., Metabacillus halosaccharovorans, V. maritimus, V. parahaemolyticus, A. flava, and Priestia aryabhattai) could degrade other carbon sources, e.g., mannose, lactose, galactose, maltodextrin, or sucrose. Our modeling showed a well-represented distribution of short-chain fatty acid production, i.e., most isolates fermented sugars into short-chain fatty acids and organic acids, including acetate, succinate, and ethanol (Additional file 1: Fig. S21).

Bacterial metabolites modulate L. marina development and longevity

Focusing on HQbiome-1 (“beneficial” to L. marina: 29 isolates from HQbiome) and HQbiome-2 (“detrimental” to L. marina: 25 isolates from HQbiome), we identified bacterial metabolic pathways associated with L. marina development through random forest regression analysis (Fig. 4a; Additional file 4: Table S3E). Especially, several bacterial metabolic pathways were significantly positively correlated with L. marina development, including ubiquinol (Coenzyme Q (CoQ): CoQ8, CoQ9, CoQ7, and CoQ10 (PWY-6708, PWY-5856, PWY-5855, and PWY-5857)), heme b (HEME-BIOSYNTHESIS-II) and acetaldehyde biosynthesis pathways (PWY66-21 and PWY-6333 with key metabolite: acetyl-CoA and acetaldehyde), L-cysteine degradation pathway (PWY-5329 with key metabolite: pyruvate), NAD phosphorylation pathways (NADPHOS-DEPHOS-PWY and NADPHOS-DEPHOS-PWY-1 with key metabolite: NAD+, NADP+, and NADPH), and N-end rule pathways (PWY-7801 and PWY-7802 with key metabolite: N-terminal arginyl-protein). Besides, several bacterial lipids metabolic pathways, including oleate and α-Kdo-(2- > 4)-α-Kdo-(2- > 6)-lipid IVA (PWY-7664 and PWY-8074 with key metabolite: oleate and α-Kdo-(2- > 4)-α-Kdo-(2- > 6)-lipid IVA), and some cofactors biosynthesis pathways, e.g., tetrahydromonapterin (PWY0-1433) and molybdenum cofactor biosynthesis (PWY-8171), were also significantly positively associated with L. marina development. The isolates belonging to Actinobacteria, Bacteroidetes, and Firmicutes lack these pathways, which might cause L. marina developmental delay. In addition, we found that five bacterial metabolic pathways were negatively correlated with L. marina development, notably, biosynthesis of pyridoxal 5’-phosphate (PLP, PWY-6466), which was reported to be critical for bacterial pathogenicity [34]. The other four negatively correlated metabolic pathways were PWY0-1517, PWY-7247, PWY-6383, and PWY-5532, respectively (Additional file 4: Table S3E).

[IMAGE OMITTED: SEE PDF]

Using indicator analysis, we identified metabolic pathways that were significantly enriched in the HQbiome-1 (Cluster1) and the HQbiome-2 (Cluster2) (Additional file 1: Fig. S22a; Additional file 4: Table S3D). Of the 65 pathways that were significantly enriched in HQbiome-1, prominent enriched functional pathways included cofactor metabolism (26.15%), lipid metabolism (21.54%), amino acid metabolism (9.23%), and nucleotide metabolism (9.23%) (Additional file 1: Fig. S22b). Pathways enriched in HQbiome-1 support the biosynthesis of lipids, including cyclopropane fatty acid (PWY0-541), (5Z)-dodecenoate (PWY0-862 and PWY-7858), oleate (PWY-7664), 10,13-epoxy-11-methyl-octadecadienoate (PWY-7691), palmitoleate (PWY-6282), and phosphatidylcholine (PWY-6826). The 61 significantly enriched pathways in HQbiome-2 were largely associated with cofactor metabolism (40.98%), carbohydrates metabolism (11.48%), secondary metabolite metabolism (9.84%), and energy metabolism (8.20%). Secondary metabolite metabolism pathways enriched in HQbiome-2 included polyketide biosynthesis (PWY-12), sugar alcohol degradation (MANNIDEG-PWY), S-containing secondary compound biosynthesis (PWY-8001), N-containing secondary compound biosynthesis (PWY-6852), and triterpenoid biosynthesis (PWY-6095 and PWY-6098).

To ask if any of the key metabolites identified from the above random forest analysis could modulate L. marina development, we evaluated the single metabolite effect through dietary supplementation. We found that dietary supplementation of heme and acetyl-CoA to Alkalihalobacillus algicola HQB-10, significantly promoted the development of L. marina in day 4 to day 6; and day 5 to day 9, respectively (Fig. 4c; the details of statistical analysis for each day were represent in Additional file 4: Table S3N). Although we observed accelerated development trend in L. marina fed HQB-10 with dietary supplementation of CoQ, which did not yield significant statistics (Fig. 4c). In addition, we found that dietary supplementation of heme, acetyl-CoA, and acetaldehyde to Metabacillus halosaccharovorans HQB-138, significantly accelerated L. marina development in days 4, 5, 7, to 10, and day 5, respectively (Fig. 4c; the details of statistical analysis for each day are represented in Additional file 4: Table S3N). Unexpectedly, we observed that dietary supplementation of oleate to both HQB-10 and HQB-138, significantly caused L. marina developmental delay (Fig. 4c, d). By contrast, dietary supplementation (e.g., CoQ, heme, acetyl-CoA, acetaldehyde, and pyruvate) did not significantly affect L. marina development fed Paenisporosarcina quisquiliarum HQB-263 or Staphylococcus warneri HQB-123, both of which lack the corresponding metabolic pathway. Unexpectedly, we found that oleate supplementation significantly attenuated L. marina growth rate or led to death (Fig. 4c, d and Additional file 1: Fig. S23). Moreover, as expected, we found that dietary supplementation of pyridoxal 5’-phosphate (PLP, B6) significantly slowed the development of L. marina in day 3 and day 4 when fed OP50, as well as in day 4 when fed HQB-5 (Fig. 4e).

To ask whether any of the key metabolites biosynthesis gene identified from the above random forest analysis, could modulate L. marina development, we generated 2 knockout mutants for each of the three genes (ubiG, hemF, and adhE), which plays pivotal roles in CoQ (PWY-6708), heme b (HEME-BIOSYNTHESIS-II), and acetyl-CoA (PWY66-21) biosynthesis, respectively, in E. coli OP50 background using homologous recombination methods [35]. Compared to wild-type animals, we found that both mutant strains of ΔubiG, ΔhemF, and ΔadhE (2 mutant strains for each gene) significantly attenuated L. marina development in day 3, respectively (Additional file 1: Fig. S24 and Fig. S25). Notably, mutations in ΔubiG significantly cased L. marina developmental delay in day 3 to 6, and day 3 to 10, respectively (Additional file 1: Fig. S25). These data suggest that ubiG, hemF, and adhE play critical roles in biosynthesis of CoQ, heme b, and acetyl-CoA to promote the development of L. marina. However, dietary supplementation of the corresponding CoQ, heme, and acetyl-CoA could not rescue the slow growth of ΔubiG, ΔhemF, and ΔadhE mutants (Additional file 1: Fig. S24 and Fig. S25), suggesting cryptic genetic background might contribute to L. marina development together with dietary supplementation of single bacteria metabolite.

We additionally characterized the metabolic interactions between all 73 HQbiome bacterial isolates and their respective effects on L. marina developmental rate (Additional file 1: Fig. S26; Additional file 4: Table S3F). Notably, most metabolic pathways associated with L. marina developmental rate (Fig. 4a) were also captured with analysis using whole HQbiome isolates (Additional file 1: Fig. S26), while 7 new pathways were identified in this new analysis, including NAD salvage, L-histidine degradation I, octane oxidation, NAD salvage pathway III (to nicotinamide riboside), thiosulfate disproportionation IV (rhodanese), 5'-deoxyadenosine degradation I, and L-tryptophan degradation IV (via indole-3-lactate).

In addition, we found that L. marina lifespan was significantly negatively correlated with glycerol degradation I pathway (PWY-4261), while positively correlated with dimethylsulfoniopropionate degradation I pathway (PWY-6046) (Fig. 4b; Additional file 4: Table S3K; the results of lifespan assays are shown in Additional file 1: Fig. S27). Notably, as expected, we found that 1% glycerol supplementation significantly extended the lifespan of L. marina fed S. algae HQB-5 but not OP50 (Fig. 4f, Additional file 1: Fig. S28). Through indicator analysis, we assessed the metabolic pathways that were significantly correlated with the longevity of Shewanella spp. and Vibrio spp. We found that 27 and 68 pathways were significantly enriched in Shewanella or Vibrio respectively (Additional file 1: Fig. S22c; Additional file 4: Table S3H). Of the 27 pathways that were significantly enriched in Shewanella, significant ones included energy metabolism (33.33%), cofactor metabolism (18.52%), and sulfur metabolism (11.11%). While 68 pathways enriched in Vibrio mainly contained cofactor metabolism (22.06%), amino acid metabolism (20.59%), carbohydrate metabolism (14.71%), and energy metabolism (14.71%) (Additional file 1: Fig. S22d).

Bacterial metabolic pathways associated with C. elegans development and longevity

Through random forest regression analysis, we identified bacterial metabolic pathways including GDP-mannose biosynthesis (PWY-5659), aerobic respiration III (PWY-4302), heme b biosynthesis I (HEME-BIOSYNTHESIS-II), and L-cysteine degradation III (PWY-5329) were significantly correlated with C. elegans development (Fig. 5a; Additional file 4: Table S3G), which were the only four pathways that also significantly correlated with L. marina growth (Fig. 4a; Additional file 4: Table S3E, G). Notably, most bacterial metabolic pathways associated with animal development were largely different between L. marina and C. elegans (Fig. 4a and Fig. 5a). Especially, we found that several bacterial metabolic pathways were significantly positively correlated with C. elegans development, including GDP-mannose biosynthesis (PWY-5659 with key metabolite: GDP-α-D-mannose), superoxide radical degradation (DETOX1-PWY), three amino acid degradation pathways (β-alanine: PWY-1781; L-alanine: ALADEG-PWY; and L-cysteine: PWY-5329 with key metabolites: pyruvate and acetyl-CoA), aerobic respiration (PWY-4302 with key metabolite: ubiquinol), heme b biosynthesis pathways (HEME-BIOSYNTHESIS-II and HEMESYN2-PWY with key metabolites: heme b), and branched-chain fatty acid biosynthesis pathways (PWY-8173, PWY-8174, and PWY-8175 with key metabolites: 12-methyltetradecanoate, 13-methyltetradecanoate, 14-methylhexadecanoate, and 15-methylhexadecanoate) (Fig. 5a; Additional file 4: Table S3G). In contrast, UDP-N-acetyl-D-galactosamine biosynthesis II, N-acetylglucosamine degradation I, and mono-trans, poly-cis decaprenyl phosphate biosynthesis pathways were negatively correlated with C. elegans egg-laying time (Additional file 4: Table S3G). Although the above key metabolites correlated with C. elegans development were distributed to different pathways compared with L. marina (Fig. 4a; Additional file 4: Table S3E), certain metabolites such as CoQ, acetyl-CoA, and pyruvate were shared between both nematodes to promote development (Additional file 4: Table S3G). Through single metabolite dietary supplementation, we found that supplementation with acetyl-CoA significantly promoted C. elegans development fed Alkalihalobacillus hwajinpoensis HQB-118, Staphylococcus warneri HQB-123, and Paenisporosarcina quisquiliarum HQB-263 (Fig. 5c–e), and pyruvate significantly promoted C. elegans development fed HQB-118 (Fig. 5c). Using the ΔubiG, ΔhemF, and ΔadhE mutants of E. coli OP50, we observed a significant developmental delay in C. elegans when fed both mutant strains of either ΔubiG or ΔhemF (Additional file 1: Fig. S29).

[IMAGE OMITTED: SEE PDF]

Next, we examined the effects of Shewanella spp. and Vibrio spp. on C. elegans lifespan (Additional file 1: Fig. S30; Additional file 4: Table S3I). Interestingly, the lifespan of C. elegans fed Shewanella spp. and Vibrio spp. was not all significantly decreased compared with E. coli OP50 (Additional file 1: Fig. S30c; Additional file 4: Table S3J). The C. elegans fed strain HQB-136 belonging to Shewanella spp. showed the longest lifespan, while this strain did not support C. elegans growth (Additional file 1: Fig. S30a, b). Additionally, we found that C. elegans lifespan was negatively correlated with adenine salvage (PWY-6610) (Fig. 5b, Additional file 4: Table S3L).

Gut bacterial colonization and bacterial impacts on L. marina and C. elegans physiology

To characterize the gut microbiota and its potential effects on animal development and lifespan, we generated germ-free L. marina inbred line worms by bleaching eggs to obtain synchronized first larval stage (L1), which were then fed the HQbiome bacterial mixture (proportional mixture of each 73 strains; Additional file 1: Fig. S31a, Additional file 5: Table S4A) on SW-NGM plates. We found that stable L. marina gut colonization was observed on day 3 of adulthood (160 h post L1; Fig. 6a), and the gut microbiome was dominated by Shewanella (especially S. algae) and Vibrio on day 3 and day 5 of adulthood, resembling the surrounding bacterial lawn (Fig. 6a; Additional file 5: Table S4B). We next used differential abundance analysis to assess variations at genus level across different time points (Additional file 1: Fig. S32). There was a notable enrichment in the presence of Shewanella and Galactobacter when L. marina grown on HQbiome on day 5 compared to day 1 (Additional file 1: Fig. S32a). When comparing the gut microbiota of L. marina to the lawn, we observed an enrichment of Aliivibrio (2.17%), Lactococcus (0.20%), and Paracoccus (0.48%) on day 3, while Parasphingorhabdus (0.22%), Bacillus (0.22%), Lactococcus (0.10%), and Paracoccus (0.21%) were enriched on day 5 compared to the other days, suggesting of selection for certain bacterial genus with low abundance (Additional file 1: Fig. S32b).

[IMAGE OMITTED: SEE PDF]

We next asked whether the mixture of the “beneficial” subset HQbiome-1 (Additional file 1: Fig. S31b; Additional file 5: Table S4A) had similar colonization characteristics compared to the HQbiome mixture. We characterized the gut microbiota of L. marina fed HQbiome-1 mixture (proportional mixture of 29 “beneficial” of 73 strains from HQbiome) and found that the gut microbiota stabilization time was similar to that of HQbiome (day 3 of adulthood; Fig. 6a, b). Although there were five Shewanella spp. in HQbiome-1 (without S. algae), strikingly, we found that L. marina gut microbiota was dominated by Vibrio (especially V. parahaemolyticus) fed HQbiome-1, resembling the bacterial lawn (Fig. 6b; Additional file 5: Table S4C). Based on differential abundance analysis, when L. marina was cultivated on HQbiome-1, Neptunomonas showed exclusive enrichment on day 3 compared to day 1, and Neptunomonas, Vibrio, and Shewanella were significantly enriched on day 5 (Additional file 1: Fig. S32e). No differentially abundant taxa were found when comparing the microbiota of L. marina and to the lawn on day 3, while Paracoccus (0.41%) and Enterovibrio (0.17%) were significantly enriched on day 5 compared to the other days, suggesting of selection for certain bacterial genus with low abundance (Additional file 1: Fig. S32f).

Given that the gut microbiota was dominated by S. algae when L. marina was fed HQbiome, and S. algae is not available in HQbiome-1, to ask whether S. algae play a major role in gut microbiota colonization, we added S. algae HQB-5 to HQbiome-1 mixture (termed HQbiome-3; Additional file 1: Fig. S31c) and characterized the bacterial composition of lawn and gut microbiota composition (Fig. 6c). Of note, the gut microbiota reverted to be dominated by Shewanella when L. marina was fed HQbiome-3 (Fig. 6c), compared with being dominated by Vibrio when fed HQbiome-1 (Fig. 6b). These data suggested that S. algae played a major role in L. marina gut microbiota colonization.

To ask if the gut microbiota modulates the physiology of L. marina, we examined the developmental rate and lifespan when L. marina grown on HQbiome (73 isolates), HQbiome-1 (29 “beneficial” isolates), HQbiome-2 (19 “detrimental” isolates), and HQbiome-3, respectively. We observed that the developmental rate of L. marina fed HQbiome-1 was significantly faster than that of HQbiome and HQbiome-2, and HQbiome significantly promoted L. marina development compared to HQbiome-2, while the developmental rate was similar when fed HQbiome-1 and HQbiome-3 (Fig. 6d). Interestingly, we found that the lifespan of L. marina fed HQbiome was significantly longer than that of HQbiome-1 and HQbiome-2, while HQbiome-1 and HQbiome-2 led to similar lifespan (Fig. 6e, f). Of note, we found that HQbiome-3 significantly shortened L. marina mean lifespan compared to HQbiome-1 feeding (Fig. 6e), indicating that S. algae played a critical role in L. marina longevity regulation.

To test whether there are similar gut microbiota characterization and bacterial effects on the development and lifespan between L. marina and its terrestrial relative C. elegans, we grew C. elegans on HQbiome, HQbiome-1, HQbiome-2, and HQbiome-3, respectively. Similar to L. marina, we observed that C. elegans gut colonization could be stabilized on day 5 of adulthood (155 h post L1) fed either HQbiome or HQbiome-1 (Fig. 7a, b). In contrast to the observation that L. marina gut microbiota was similar in composition to their bacterial lawn (Fig. 6a–c), we found that C. elegans gut microbiome composition exhibited selectivity from their bacterial lawn (Fig. 7a–c and Additional file 1: Fig. S33). Especially, the gut microbiota of C. elegans grown on HQbiome was dominated by Shewanella and Pseudomonas on day 5 and day 7 of adulthood (Fig. 7a), while C. elegans fed HQbiome-1 was dominated by Pseudomonas and Vibrio (Fig. 7b). Notably, we observed that C. elegans gut microbiota was dominated by Shewanella and Pseudomonas when animals grown on HQbiome-3 (Fig. 7c), indicating that S. algae played a critical role in C. elegans gut microbiota colonization. Though differential abundance analysis, when C. elegans grown on HQbiome, Microbacterium, Arthrobacter, and Shewanella exhibited significant enrichment on day 5 compared to day 1, whereas no species showed significant enrichment between day 5 and day 7, suggesting the stability of the microbiota of C. elegans from day 5 (Additional file 1: Fig. S32c). In contrast to the lawn, Vibrio (1.24%) and Arthrobacter (0.28%) were significantly enriched in the gut microbiota of C. elegans on day 5, and Vibrio (0.59%) continued to be enriched on day 7 (Additional file 1: Fig. S32d). Comparing different days of the gut microbiota of C. elegans grown on HQbiome-1, we did not identify any differentially abundant taxa. Enterovibrio (19.26%) was exclusively enriched in the gut microbiota of C. elegans compared to the lawn on days 3, 5, and 7, while Pseudomonas (85.44%) and Saccharospirillum (0.25%) were significantly enriched on day 5, when C. elegans grown on HQbiome-1 (Additional file 1: Fig. S32g). However, no differentially abundant taxa were found when comparing the gut microbiota of C. elegans and to the lawn on day 3 when C. elegans grown on HQbiome-3, suggesting lack of selection. Similar to L. marina, we observed that HQbiome-1 significantly shortened C. elegans egg-laying time compared to HQbiome and HQbiome-2, whereas HQbiome-2 significantly delayed egg-laying time compared to HQbiome (Fig. 7d). Of note, C. elegans exhibited significant shorter egg-laying time fed HQbiome-3 compared to HQbiome-1, suggesting that S. algae played an essential role in promoting C. elegans development. By contrast, we found that the lifespan of C. elegans fed HQbiome-2 was significantly increased compared to HQbiome, HQbiome-1, and HQbiome-3, indicating that tradeoffs occurred between C. elegans developmental rate and lifespan (Fig. 7d–f). In addition, we simulated the growth and interactions within bacterial communities and found that the H2S concentration was twice as high in HQbiome-2 compared to HQbiome-1 (Additional file 1: Fig. S34 and Fig. S35). This increased concentration could contribute to the observed lifespan extension [36]. Collectively, these data showed that the impacts of three bacterial mixtures (HQbiome, HQbiome-1, and HQbiome-3) on C. elegans developmental rate were in similar trends to that of marine nematode L. marina, whereas the lifespan of nematodes might be not only determined by the dominant gut bacteria but also by the interactions between microbiota members in a complex gut bacterial ecosystem.

[IMAGE OMITTED: SEE PDF]

Discussion

In the present study, we applied both cultivation-independent and cultivation-based approaches to characterize the environmental sedimentary bacterial composition in Huiquan Bay, where marine nematode L. marina was collected. We demonstrated that four phyla were dominant, including Proteobacteria, Bacteroidetes, Cyanobacteria, and Actinobacteria, in accordance with the prominent habitat microbes of terrestrial nematode C. elegans (Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria) [15]. At the family level, Flavobacteriaceae was abundant in L. marina native habitats (8.86%), similar to that in C. elegans habitats (3.16%) [15]. Among the total OTUs from the cultivation-independent full-length 16S rRNA gene sequencing, only 1.43% (0.43% relative abundance in sediment samples) were recovered in our cultivation-dependent bacterial isolation, which was coincident with earlier studies that < 1% of marine environmental bacteria were culturable [37]. Notably, through the cultivation-dependent approach, we isolated 539 habitat bacteria, including 28 new species.

Both marine nematode L. marina and its terrestrial relative C. elegans belong to the same family, Rhabditidae [26]. The comparative investigation of L. marina and C. elegans could provide valuable insights into the adaptations of nematodes to the bacterial environments [27]. We examined the effects of 448 marine native bacteria isolates, which were taxonomically similar to the habitat microbiome of C. elegans at the phylum level, on the development of C. elegans. We found that approximately 82% of single isolates were able to support the C. elegans growth, while 18% of isolates significantly delayed C. elegans development (Additional file 1: Fig. S8). By contrast, we observed that approximately 52.34% of 448 native bacteria supported L. marina development, while 47.66% caused L. marina developmental delay or death (Fig. 2a). Consistent with a previous report that about 80% of habitat-bacterial isolates support C. elegans growth, our data indicated that C. elegans could survive in a wider range of bacterial community niches. While a greater proportion of bacterial isolates promote C. elegans growth, upon evaluating the impacts of all 448 individual bacteria on the developmental rate of L. marina and the egg-laying time of C. elegans, we found a positive correlation in the effects of bacteria on nematodes’ development between L. marina and C. elegans (Additional file 1: Fig. S9a, b), indicating the conservational role of marine native bacteria in both nematode’ physiology. However, this correlation is not significant within certain phylum (e.g., Betaproteobacteria; Additional file 1: Fig. S9c; Rs = 0.21, P = 0.79).

Among the 73 HQbiome isolates, 66% (40% “beneficial” and 26% “intermediate”) supported L. marina development, while 34% of the isolates were detrimental to L. marina development (Fig. 2b). Similar to the effects of 448 marine native bacteria on C. elegans development, 77% of HQbiome isolates (68% “beneficial” and 9% “intermediate”) supported C. elegans growth, while 23% of the isolates slowed C. elegans development. In line with 448 isolates’ effects on L. marina growth rate and the egg-laying time of C. elegans, HQbiome isolates affect the development of both nematodes in a similar way (Additional file 1: Fig. S11).

Based on our reconstructed metabolic models from WGS of 73 HQbiome isolates and their effects on the development of L. marina, we found that several metabolic pathways were significantly positively correlated with L. marina development rate (Fig. 4a). Essential metabolites from the above pathways included CoQ, heme b, acetyl-CoA, acetaldehyde, pyruvate, and oleate. In contrast, we found that the siderophores, pyridoxal 5’-phosphate, and terpenes were negatively correlated with animal growth (Fig. 4a and Additional file 1: Fig. S15d). Notably, the deletion of the essential genes of CoQ, heme b, and acetyl-CoA attenuated L. marina development (Additional file 1: Fig. S24 and Fig. S25), and the supplementation of heme, and acetyl-CoA significantly promoted L. marina development fed A. algicola HQB-10, and dietary supplementation of heme, CoA, and acetaldehyde significantly promoted L. marina growth fed M. halosaccharovorans HQB-138 (Fig. 4c, d), in accordance with lack of the corresponding metabolic pathways in both bacteria. Unexpectedly, we found that oleate supplementation caused L. marina developmental delay or death (Fig. 4c, d and Additional file 1: Fig. S23). CoQ plays an essential role in mitochondrial energy metabolism [38] and antioxidant activity [39, 40]. Mutations in C. elegans clk-1 led to synthetic CoQ defect, showing development delay, reduced respiration, and increased life span [41,42,43,44]. Heme is a readily bioavailable source of iron for animal life [45]. Heme b was reported to be present in the respiratory complexes to facilitate the movement of electrons through the electron transport chain, and to stabilize ubisemiquinone radicals and decrease the production of reactive oxygen species (ROS) [46]. Given heme could not be synthesized by many nematodes, such as C. elegans, which requires environmental heme for growth and development [45], our data indicated that L. marina could not synthesize heme b as well. The production of vitamin B6 was reported to be essential for bacterial pathogenicity to C. elegans [34], which might delay L. marina development due to increased susceptibility to pathogen bacteria. As expected, PLP supplementation significantly attenuated the development of L. marina fed S. algae HQB-5 (Fig. 4e), which lacks PLP biosynthesis pathways.

Only four bacterial metabolic pathways, including GDP-mannose biosynthesis, aerobic respiration III, heme b biosynthesis I, and L-cysteine degradation III, were significantly positively correlated with L. marina developmental rate and C. elegans egg-laying time. However, most of the correlated metabolic pathways were not common to both nematodes. In contrast to the pro-development role of heme in L. marina, heme dietary supplementation did not affect C. elegans egg-laying time. The mechanisms by which single bacterial metabolite regulates nematode development deserved to be further explored.

Our study demonstrated that the glycerol degradation I pathway was significantly negatively correlated with L. marina lifespan, it has been reported that the dietary supplement of glycerol extended the lifespan in both yeast and rotifers [47, 48]. We speculated that the bacteria that lack glycerol degradation competence may supplement L. marina with more glycerol and extend the nematode lifespan. In line with our hypothesis, we observed that dietary supplementation of 1% glycerol significantly promoted the longevity of L. marina fed S. algae HQB-5 but not OP50 (Fig. 4f). By contrast, 1% glycerol supplementation did not extend C. elegans lifespan fed either OP50 or S. algae HQB-5 (Additional file 1: Fig. S28a, c). Several studies in C. elegans have uncovered that glycerol supplementation shortened worm lifespan [49, 50]. The underpinning mechanisms by which glycerol supplementation regulates L. marina longevity deserved to be further explored.

A recent report showed that C. elegans gut microbiota stabilized at day 3 of adulthood [14], while our data revealed that the gut microbiota stabilization occurred at 160 and 150 h post L1 larvae for L. marina and C. elegans, respectively (i.e., at day 3 of adulthood for L. marina and day 5 of adulthood for C. elegans), suggesting that distinct environmental bacterial composition may affect the stabilization time of gut microbiota colonization. It was reported that C. elegans exhibited three distinct microbiota selection strategies (strong selection, medium selection, and without selection), and the wild-type N2 worms showed medium selection [14]. In agreement with this report, we found that gut microbiota colonization in C. elegans N2 worms exhibited medium selection (Fig. 7a, b and Additional file 1: Fig. S33). It was described that C. elegans gut microbiota selection was determined either by insulin receptor gene daf-2 or the TGFβ/BMP-like ligand gene dbl-1 [14, 23]. However, we found that the gut microbiota of L. marina was similar in composition to the bacterial lawn, suggesting a lack of selective filtering (Fig. 6a, b and Additional file 1: Fig. S33), implying that certain selection might occur in daf-2 and dbl-1 and other host genes in L. marina genome [14, 23]. The gut microbiota was dominated by S. algae when L. marina grown on HQbiome, and Vibrio was the dominant microbe when fed on HQbiome-1 where S. algae was unavailable. Strikingly, the gut microbiota reverted to be dominated by Shewanella spp. when L. marina fed HQbiome-3 (HQbiome-1 plus S. algae). These data suggested that S. algae played a major role in L. marina gut microbiota colonization. Notably, we found that S. algae also played an essential role in C. elegans gut microbiota colonization. Further studies to elucidate the casual interactions between habitat microbiome, gut microbiota, and host biology will provide novel insights into the conservation and management of ecosystems.

Of note, the development of nematodes fed different bacteria mixtures were largely positively correlated, i.e., HQbiome-2 caused both nematode developmental delays compared to all the other bacterial mixtures, although the effect of S. algae in HQbiome-3 on the development of L. marina and C. elegans are different (Fig. 6d and Fig. 7d). However, the lifespan of L. marina and C. elegans fed different bacteria mixtures were largely different, i.e., HQbiome-2 caused totally opposite effects on the life span of both nematodes, and the effect of S. algae in HQbiome-3 on the life span of L. marina and C. elegans are different (Fig. 6e and Fig. 7e). Bacteria are food sources for many nematodes such as L. marina and C. elegans, thus, environmental bacteria might provide essential nutrients that support animal host development. We found that certain bacteria such as S. algae showed different effects on animal physiology when nematodes host fed different bacterial mixtures, indicating that there are complex interactions between bacteria species. These data suggested that the same bacteria mixtures might play different roles in the lifespan between L. marina and C. elegans, suggesting different host–microbe interactions.

Nematodes have evolved through numerous habitat transitions between ocean and land followed by moderate diversification, it was described that L. marina evolved from a transition from terrestrial to marine habitats [51]. Further microbiome-mediated comparative functional analysis between marine nematode L. marina and its terrestrial relative C. elegans, will permit future microbiota manipulation experiments to decipher how habitat microbes shape animal host fitness and their evolutionary transition trajectories.

Conclusions

Hundreds of marine isolates collection and their effects on nematode physiology presented here will open avenues for further mechanistic studies about microbiome–host interactions. Given that both bacteria and marine nematodes are dominant taxa in sedimentary ecosystems, the resource presented here will provide novel insights to identify mechanisms underpinning how habitat bacteria affect nematode biology in a more natural context. Our integrative approach will also provide a novel marine microbe–nematode framework to elucidate microbiome-mediated functions and their contribution to healthy ecosystems, in the context of global climate change.

Methods

Cultivation-independent full-length 16S rRNA gene profiling

Eighty samples were collected from Huiquan Bay, Qingdao, China, from June 2020 to June 2021. Full-length 16S rRNA gene sequencing was applied to the 80 sedimentary samples. Extraction and purification of total DNA from 0.25 g of sediment samples (wet weight) were carried out with the TGuide S96 Magnetic Soil/Stool DNA Kit (Tiangen biochemical technology, China) according to the manufacturer’s protocols. PCR amplifications were performed using primers flanking full-length of the bacterial 16S rRNA genes (27F: 5′-AGRGTTTGATYNTGGCTCAG-3′ and 1492R: 5′-TASGGHTACCTTGTTASGACTT-3′). The 16S full-length reaction system consists of a 30 μL mixture, which contains 1.5 μL of genomic DNA, 10.5 μL of NFW (Hydration Volume), 15 μL of KOD One TM PCR Master Mix (Biotechnology Co., Ltd., China), and 3 μL of the barcode primer pair. The PCR parameters for the detection of the 16S full-length region primer samples were as follows: 95 °C pre-denaturation for 2 min; 30 cycles of 98 °C denaturation for 10 s, 55 °C annealing for 30 s, and 72 °C extension for 90 s; and 72 °C extension for 2 min. Then the PCR products were purified, quantified, and homogenized to construct a PacBio HiFi sequencing library using SMRT bell Express Template Prep Kit 2.0 (Pacific Biosciences, US), which was sequenced on PacBio Sequel II platform (Pacific Biosciences, US). Raw reads were processed using the SMRT Link v 8.0 software with the minPasses = 5 and minPredictedAccuracy = 0.9 to obtain the circular consensus sequencing (CCS) reads. Subsequently, CCS reads were then demultiplexed using Pacific Biosciences barcode demultiplexer (lima) to identify barcode sequences. CCS reads containing no primers and those reads beyond the length range (1200–1650 bp) were discarded through the recognition of forward and reverse primers and quality filtering using the Cutadapt quality control process (version 2.7) [52]. Quality filtering, dereplication, learning the dataset-specific error model, amplicon sequence variants (ASVs) inference, and chimera removal were performed using DADA2 (dada2 denoise-ccs) pipeline [53] in QIIME2 v2022.2.0 [54]. After quality control filtering, the high-quality merged sequences were clustered into operational taxonomy units (OTUs) using the VSEARCH algorithm (vsearch cluster-feature-de-novo) [55] with a threshold of 97% pair-wise nucleotide sequence identity. Taxonomy annotation of the OTUs was performed based on the VSEARCH-based consensus classifier [55] using the EzBioCloud database. Statistical analyses and graphical visualization were obtained using R v4.1.3. PICRUSt2 v2.5.0 [56] was used to predict functional pathways in modules of microbial taxa. Additional meta-information on predicted PICRUSt2 output was obtained from the MetaCyc pathway Database [57].

Cultivation-dependent bacterial isolation

Intertidal sediments from Huiquan Bay, Qingdao, China (36.05 N, 120.33 E) were collected every about 2 months over about a 3-year period (August 2019 to April 2022; Additional file 1: Fig. S1). Especially, about 30 g sediment samples were collected using 50-mL sterile centrifugal tubes, then 10 g sediment was transferred to a new 50-mL centrifugal tubes with 30 mL sterile seawater and was vortexed. For large-scale isolation, 100 μL of supernatant in serial dilutions were plated onto a combination of 10 different media with various culture conditions (see “Media” section). Culture plates were incubated at 28 °C overnight. A total of 3822 bacterial isolates were picked randomly for DNA amplification and sequencing of their 16S rRNA gene. The 16S rRNA gene sequences were PCR amplified using bacterial primers 27F (5′-AGRGTTYGATYMTGGCTCAG-3′) and 1492R (5′-RGYTACCTTGTTACGACTT-3′), and after PCR-product cleaning, Sanger sequencing was performed with primer 27F. EzBioCloud and NCBI databases were applied to search for the most closely related species. For storage, cryo-stocks were prepared by mixing freshly grown cultures 1:1 in 30% (v/v) glycerol and stored at − 80 °C.

Media

Ten different media were modified from several base-type that favored the growth of different bacteria:

1. (1)

2216E: Tryptone 5 g, Yeast extract 1 g, agar 20 g, seawater 1 L.

2. (2)

2216E+: Tryptone 5 g, Yeast extract 1 g, NaCl 30 g, agar 20 g, seawater 1 L.

3. (3)

1/10 2216E: Tryptone 0.5 g, Yeast extract 0.1 g, agar 20 g, seawater 1 L.

4. (4)

1/10 2216E+: Tryptone 0.5 g, Yeast extract 0.1 g, NaCl 30 g, agar 20 g, seawater 1 L.

5. (5)

R2A: Tryptone 0.25 g, Peptone 0.25 g, Casein Hydrates 0.5 g, Yeast extract 0.5 g, Soluble starch 0.5 g, K2HPO4 0.3 g, MgSO4 0.1 g, Glucose 0.5 g, Sodium pyruvate 0.3 g, agar 12 g, ddH2O 1 L, pH 7.2 ± 0.2.

6. (6)

R2A+: Tryptone 0.25 g, Peptone 0.25 g, Casein Hydrates 0.5 g, Yeast extract 0.5 g, Soluble starch 0.5 g, K2HPO4 0.3 g, MgSO4 0.1 g, Glucose 0.5 g, Sodium pyruvate 0.3 g, agar 12 g, NaCl 30 g, ddH2O 1 L, pH 7.2 ± 0.2.

7. (7)

1/3 R2A: Tryptone 0.08 g, Peptone 0.08 g, Casein Hydrates 0.17 g, Yeast extract 0.17 g, Soluble starch 0.17 g, K2HPO4 0.1 g, MgSO4 0.03 g, Glucose 0.17 g, Sodium pyruvate 0.1 g, agar 12 g, ddH2O 1 L, pH 7.2 ± 0.2.

8. (8)

1/3 R2A+: Tryptone 0.08 g, Peptone 0.08 g, Casein Hydrates 0.17 g, Yeast extract 0.17 g, Soluble starch 0.17 g, K2HPO4 0.1 g, MgSO4 0.03 g, Glucose 0.17 g, Sodium pyruvate 0.1 g, agar 12 g, NaCl 30 g, ddH2O 1 L, pH 7.2 ± 0.2.

9. (9)

SCA: Soluble starch 10 g, Casein 0.3 g, K2NO3 2.0 g, NaCl 30.0 g, KH2PO4 2.0 g, CaCO3 0.02 g, MgSO47H2O 0.05 g, FeSO47H2O 0.01 g, agar 15 g, seawater 1 L.

10. (10)

LB: Tryptone 10 g, Yeast extract 5 g, NaCl 10 g, agar 15 g, ddH2O 1 L.

Nematodes maintenance

Litoditis marina was a 23rd generation inbred line generated by consecutive full-sibling crosses in our lab [26]. The Caenorhabditis elegans N2 strain utilized in this study was obtained from the Caenorhabditis Genetics Center (CGC, https://cgc.umn.edu/). C. elegans strains grown on nematode growth media (NGM), while L. marina strains were maintained on seawater nematode growth media (SW-NGM) [26]. For regular maintenance, both NGM and SW-NGM were seeded with E. coli OP50 as a food source and worms were cultured at 20 °C.

Developmental rate assay for L. marina

To assess the influence of different habitat bacterial isolates on the development of L. marina, we examined the growth rate of L. marina on single and mixed bacterial lawns. Prior to the phenotypic assays, all isolates were incubated at 28 °C until single colonies were visible. The single colonies were transferred to 30-mL culture media and incubated at 28 °C overnight. Fifteen microliters of the bacterial culture (OD600 = 0.5) were spotted into 3.5-cm SW-NGM plates. The standard food E. coli OP50 was assayed for comparison [15]. Approximately 70 synchronized L1 worms were transferred onto the plates, and the number of L4 worms was scored from day 3 to day 10. Four hundred forty eight isolates from cultivation-dependent collections were screened individually in at least one replicate. As for 73 HQbiome isolates at least three replicates were scored. Since the proportion of L4 of L. marina after 5 days post-hatching was the most relevant influencing factor in the PCoA analysis (Additional file 1: Fig. S12c, Rs = 0.9546475, P = 2.2 × 10−16), meaning that growth phenotypes could be best distinguished, we chose the proportion of L4 after 5 days post-hatching as the comparison index for further studies. The clustering of the impact of culture collection on L. marina was analyzed by k-means clustering [58] of the percentage of L4 larvae after 5 days post-hatching (k = 3, n = 1000).

Egg-laying time assay for C. elegans

To assess the influence of 448 bacteria isolates on C. elegans development, we examined the egg-laying time of wild type C. elegans N2. Approximately 35 synchronized L1 larvae were transferred to 3.5-cm-diameter NGM plates containing different isolates, and the egg-laying time was recorded from 52 to 168 h (egg-laying time > 168 h was recorded as 168 h). At least three replicates were recorded for each isolate and mixture. The clustering of the impact of culture collections on C. elegans was analyzed by k-means clustering [58] of the egg-laying time of C. elegans, respectively (k = 3, n = 1000).

Lifespan assays

The lifespan assay was performed starting at day 1 of adulthood as described previously [59]. Briefly, 3.5-cm-diameter assay plates seeded with 15 µL single or mixture bacteria were prepared every day. Forty C. elegans L4 larvae or L. marina L4 females were transferred to 3.5-cm NGM and SW-NGM, respectively. The number of live, dead, and missing worms was scored every 24 h. Worms were scored as dead if no response was detected after prodding with a platinum wire. Dead worms on the wall of the plate were scored as missing. Alive worms were transferred to fresh lawn-seeded agar plates daily. Three replicates were analyzed for each bacterial lawn or different mixtures.

Whole genome sequencing and analyses

We generated WGS for 72 bacterial isolates from HQbiome, of which, 15 were sequenced using both Oxford Nanopore Technologies (ONT) and Illumina technologies, one with both Pacific Biosciences (PacBio) and Illumina technologies. The remaining 56 strains were sequenced with Illumina only, and the standard food E. coli OP50 was also sequenced with Illumina in this study. All genomic DNA was extracted with the Sodium Dodecyl Sulfate (SDS) method [60]. For Illumina, sequencing libraries were generated using NEBNext® Ultra™ DNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer’s instructions, then sequencing was performed on the Illumina NovaSeq PE150. Low-quality reads and adaptors were trimmed with Trimmomatic v0.36 [61]. Clean data was assembled using SPAdes v3.13.0 [62], and contigs less than 1000 bp were removed. For ONT, sequencing libraries were constructed with a 10-kb insert size using the EXP-NBD104 barcoding kit (Oxford Nanopore Technologies, UK) and SQKLSK109 connection kit (Oxford Nanopore Technologies, UK). Libraries were sequenced on a Nanopore PromethION platform (Oxford Nanopore Technologies, UK). ONT data was assembled with Illumina data using SPAdes v3.13 [62], Miniasm v0.3 [63], Racon v1.4.3 [64], and Unicycler v0.4.8 [65]. For PacBio, the SMRTbell™ template library was constructed according to the manufacturer’s manual by using SMRTbell Express Template Prep Kit 2.0 (Pacific Biosciences, US) with an insert size of 10 kb, and sequenced with a PacBio Sequel platform (Pacific Biosciences, US). The complete genome assembly of the PacBio reads was performed using SMRT Link v5.0.1 software [66]. The preliminary assembly was corrected by the variant Caller module of the SMRT Link software v5.0.1 [66] using Illumina reads aligned by BWA [67]. We sequenced and assembled the E. coli OP50 genome and got a high-quality genome (8 contigs, 100% completeness, Additional file 3: Table S2F). Detailed assembly statistics for each sequenced isolate can be found in Additional file 3: Table S2F. The quality of all genome assemblies was assessed with BUSCO v5.2.2 [68] and QUAST v5.0.2 [69] (Additional file 1: Fig. S36). Of note, there is no difference between BUSCO and MEMOTE parameters (see below) between complete genomes and incomplete genomes (Additional file 1: Fig. S37). Genomes were compared using BRIG [70] for members of the same genus/family (Additional file 1: Fig. S13). Genomes were then annotated using the PROKKA v1.14.6 [71]. Pathway analysis was conducted by Eggnog-mapper v2.1.6 [72]. Average nucleotide identities were visualized using FastANI v1.33 [73] and ANIclustermap for species validation. The whole genome-based tree was prepared using ezTree [74]. To visualize the phylogenetic tree topology, the raw files were processed using iTOL v6.5.7 [75].

Phylogenetic signal analyses

The presence of a phylogenetic signal in impacts on nematode growth was tested by calculating Pagel’s λ using the package phytools [76] and by calculating Abouheif’s Cmean using the package adephylo [77]. To test whether the uneven taxonomic distribution affects the stability of this analysis, the metrics on randomly subsampled all families with more than five isolates were additionally analyzed (Additional file 1: Fig. S7). One thousand iterations for each were performed.

Metabolic network reconstructions

The 74 whole genome sequences (73 HQbiome isolates and OP50 as a control) were used to reconstruct genome-scale metabolic models using the gapseq v1.2 analysis pipeline [78]. Pathways were predicted from the MetaCyc pathway database v24.5 [57] and transporters were predicted based on the Transporter Classification Database (TCDB) [79]. In more detail, pathways and transporters were predicted by gapseq find (bitscore > 150 and a coverage of at least 75%), and the draft network was created by gapseq draft (-u 200, -l 100). Potential gaps in the network were resolved using gap-filling (-b 100). All other gapseq settings were kept at their default values. We used 2216E for marine bacteria and LB for OP50 as the growth medium in gap-filling step to allow in silico growth using flux balance analysis (FBA) [80]. Especially, the components of the LB medium can be found in the gapseq repository; the components of 2216E were inferred from https://www.dsmz.de/microorganisms/medium/pdf/DSMZ_Medium514.pdf. The specific medium composition is encoded in the models and the resulting models were used for further analysis. The medium composition and metabolic models are available in the GitHub repository. To test the quality of our models, we used the MEMOTE v0.12.0 analysis pipeline [32]. In brief, the average model quality scores for all models resulted in 77.74%. Most importantly, the key requirements for the models all reached at least 99%.

BIOLOG experiments

We further used BIOLOG GEN III plates (Biolog, Hayward, CA, US) to examine the metabolic competences of OP50 and four Vibrio strains (HQB-20, HQB-21, HQB-114, and HQB-171) following the manufacturer’s instructions. In brief, overnight cultures were centrifuged at 9000 g for 3 min to remove the growth medium. This step was followed by washing three times using Biolog’s inoculation fluid and resuspension in inoculation fluid (IF-A for OP50, IF-B for Vibrio strains) to a starting optical density (OD600) of 0.05. Inoculums were added into the BIOLOG GEN III plate wells and cultivated at 28 °C for 48 h. Absorbance values at 590 nm were measured, and substrate reduction was inferred from fold-change in absorbance [22]. Fold-changes in inoculation fluid were subtracted as background.

Carbon sources utilization modeling

The carbon source utilization by the HQbiome strains was predicted based on electron transfer from the source to electron carriers (i.e., ubiquinol, menaquinol, or NADH), which is analogous to the experimental carbon source test of the BIOLOG test described above [78]. In brief, temporary reactions to recycle reduced electron carriers as utilization indicators of the focal potential carbon source were added to each model, and the flux through these recycling reactions was set to max. If at least one of the recycle reactions showed a positive flux, the utilization of the carbon source was confirmed. An average of 73.20% of the experimental data were consistent with carbon utilization predictions, which was greater than the previous study [22, 81]. We then used the ProTraits database (confidence threshold of ≥ 0.95, http://protraits.irb.hr/data.html) to assess the carbon source metabolic capability of the strains present in the database. A 90% correct prediction (37/41) highlights the quality of the HQbiome models.

Fermentation product simulation

Carbon source utilization and metabolic by-products were predicted using the published script [78]. In detail, by-products from anaerobic metabolism were simulated using FBA coupled with minimization of total flux. Meanwhile, the maximum fermentation product release of individual metabolites across all possible FBA solutions was predicted using FVA [82]. A model was considered to be able to secrete a fermentation product if the exchange flux (i.e., outflow) through the fermentation product exchange reaction was positive [78].

Simulation of bacterial communities in silico growth

We performed dynamic and agent-based simulations of bacterial population growth and metabolic fluxes using the R package “BacArena” [83]. In summary, all metabolic models were simulated within microbial communities using the “BacArena,” which enables the dynamic simulation of individual-based metabolic models, which are separately optimized within a common growth environment. The size of arena was 40 × 40 grid cells. The gap-filling medium was employed to ascertain the initial substance concentrations within the arena. The community underwent simulation for a span of seven-time steps, and the analysis of metabolite uptake and production rates was conducted following the fourth time step.

Phylogenetic correlation and clustering of metabolic pathways

We assessed the correlation between phylogenetic relationships and metabolic pathway similarity, using pairwise comparisons of HQbiome strains [22]. To calculate phylogenetic similarities, pairwise similarities in 16S rRNA sequences were scored as percent identity using the full-length 16S rRNA gene sequences extracted from whole genome sequencing data using RNAmmer v1.2 [84]. To determine metabolic pathway similarity between isolates, metabolic networks from gapseq [78] were treated as vectors and clustered horizontally, followed by computation of Euclidean distances between vectors. Cluster similarity was estimated by average linkage and assessed via multi-scale bootstrapping (1000 replications) using pvclust [85].

Analyses of functional and phenotypic diversity between sequenced isolates

Analyses of functional diversity between sequenced isolates were conducted using metabolic networks as described above. Phenotypic diversity of the two nematodes was conducted by generating a profile of the growth phenotype of L. marina (proportion of L4 from day 3 to day 10) and C. elegans (egg-laying time), respectively. Subsequently, a distance measure based on the Euclidean distances of each pair of phyletic patterns was calculated, which allowed us to embed each genome as a data point in a metric space. PCoA was performed on this space of functional and phenotypic distances using custom scripts written in R. Pairwise functional distances within each phylum level (Proteobacteria are shown at the class level) cluster was performed by calculating the average distance between all pairs of genomes belonging to each cluster. Only phyla with at least five members are calculated.

Adaptive strategies

We classified HQbiome isolations using published universal adaptive strategy theory (UAST) criteria based on scores inferred from metabolic model data [22]. Briefly, UAST would divide heterotrophic bacteria into three broadly defined groups that share similar functional capabilities and ecological strategies: (i) “competitor” taxa that outcompete other microbial taxa for space or resources, (ii) “stress-tolerant” taxa that can persist under low-resource conditions, and (iii) “ruderal” taxa that can grow rapidly and exploit unoccupied niches. The criteria of “competitor” taxa were large genome size, antibiotics production (“Antibiotic-Biosynthesis” MetaCyc category based on metabolic models), siderophore biosynthesis (MetaCyc: “Siderophores-Biosynthesis”), and high catabolic diversity (MetaCyc: “Energy-Metabolism”). The standards of “stress-tolerant” taxa were auxotrophies (predicted by KBase [86]), including few rRNA copies, slow growth rates in 2216E, and exopolysaccharides production (MetaCyc: PWY-6773, PWY-6655, PWY-6658, PWY-1001, PWY-6068, PWY-6082, and PWY-6073). The characteristics of “ruderal” taxa included low catabolic diversity (MetaCyc: “Energy-Metabolism”), multiple rRNA copies, and a high growth rate in 2216E. The strategies of each isolate were assessed based on the relative score of each classification using quantile value. An isolate was assumed to follow the strategy, for which it produced the highest score. A mixed strategy was assumed if two strategies had the same score.

Metabolic pathways enrichment analysis

Detection of enriched metabolic pathways between different groups was carried out using a custom R script. In brief, the pathway table from gapseq (constant columns removed) was used as input for indicator species analysis implemented as the multipatt() function in the R indicspecies package [87]. The statistical significance of enriched pathways frequency was estimated by permutation using 9999 permuted groupings of the genomes. P values were corrected for multiple testing using FDR with a value of FDR < 0.05 being considered significant. Pathways enriched in HQbiome-1 (cluster1) and HQbiome-2 (cluster2) were displayed as a heatmap using the superheat package in R. Functional category assignments for pathways were derived from MetaCyc hierarchies and manually curated, and quantified in R.

Random forest regression analysis

We analyzed the association between phenotypic measurements (i.e., the developmental rate and lifespan) and bacterial metabolic capacities using Spearman rank correlation and random forest regression analysis as previously described [22]. A random forest regression model was performed to select features via a permutation-based score of importance using the R package VSURF v1.1.0 [88] by default settings (ntree = 2000, mtry = p/3). The predicted metabolic pathways from gapseq were employed as predictors. The mean growth rate of L. marina on day 5, or the egg-laying time of C. elegans, based on three or more replicates, served as the response variables for regression.

Single bacterial metabolite supplementation

We sought to validate the impact of pivotal key metabolites identified from the random forest model on L. marina and C. elegans development and longevity. The stock solutions of heme (1 mM), acetyl-CoA (10 mM), acetaldehyde (50 mM), pyruvate (25 mM), and glycerol (50%) were prepared in ddH2O and filter-sterilized. Stock solutions of coenzyme Q10 (100 μM) and oleate (5 mM) were prepared in DMSO and filter-sterilized. A 1 mM stock solution of vitamin B6 was made in M9 and filter-sterilized. The stock solution was diluted to final concentrations as follows: 10 μM coenzyme Q10, 20 μM heme, 0.25 mM acetyl-CoA, 5 mM acetaldehyde, 2.5 mM pyruvate, 500 μM oleate, 1 μM vitamin B6, and 1% glycerol. Three concentrations were tested to avoid the toxicity (Additional file 1: Fig. S23 and Additional file 1: Fig. S38).

Construction of the mutant strains

Deletion of ubiG, hemF, and adhE from the chromosome of WT was performed using homologous recombination methods, based on the λ red recombinase system [35]. The plasmids pKD3, pKD46, and pCP20, which were utilized in mutant construction, were generously provided by the laboratories of Guang Zhao at Shandong University and Min Xu at the Tianjin Institute of Industrial Biotechnology, Chinese Academy of Sciences. Briefly, primers containing a 50-bp region homologous to the 5’ and 3’ extremities of the target genes were used to amplify the chloramphenicol (cat) cassette from plasmid pKD3 by PCR amplification (Additional file 1: Table S6). Plasmid pKD46-containing E. coli OP50 was grown in ampicillin-containing LB broth at 30 °C until OD600= 0.2. Then, 30 mM L-arabinose (BBI, Shanghai, China) was used to induce expression of the λ Red proteins from pKD46 until OD600 = 0.6. Electrocompetent cells were prepared by concentrating 100-fold and washing 3 times with ice-cold, sterilized 10% glycerol. The purified cat cassettes were electroporated using Bio-Rad Micropulser (Bio-Rad, Inc., Hercules, CA, USA) to generate homologous recombination. Recombinant bacteria were screened, and cat-resistant colonies were grown on cat-containing LB agar at 42 °C overnight. After incubation, one colony per clone was selected and subcultured in cat-containing LB agar at 42 °C. Allelic replacements of target genes by the cat cassette were verified by PCR screening using specific primers (Additional file 1: Table S6) and DNA sequencing. Finally, the transductants were transformed with plasmid pCP20 to eliminate the cat cassette by FLP-mediated excision, followed by growth of transformants at 42 °C for plasmid curing. The resulting E. coli OP50 mutant strains were verified by PCR, and stored at −80 °C.

Preparation of microbiome mixtures

In brief, all HQbiome strains were grown on 2216E or R2A plates and incubated at 28 °C until single colonies were visible. The colonies on the plate were then grown overnight at 28 °C and 220 rpm shaking in 15-mL centrifuge tubes (ExCell Bio, China) filled with 10 mL 2216E/R2A. Cultures were harvested by centrifugation at 4000 g for 10 min, adjusted to an OD600 of 0.1 in sterile filtered M9 buffer. Inoculums of bacterial mixtures were prepared by combining equal volumes of each isolate, which was then seeded (15 μL) on NGM and SW-NGM plates. Seeded plates were grown overnight at room temperature before use. Gut microbiome samples were collected following the previously published protocol [14].

Gut microbiome sample collections

Gut microbiome samples were collected following the previously published protocol [14]. Around 150 germ-free and synchronized L1 worms were transferred in triplicate onto the HQbiome lawn at 20 °C. Nematode populations and bacterial lawns were harvested at days 1, 3, 5, and 7 of adulthood. On the sampling day, worms were picked from the lawn using a sterile platinum wire to a sterilized 1.5-mL centrifuge tube (Biosharp, China) filled with 1 mL M9 buffer (M9 + 0.01% Triton X-100). The tubes were centrifuged at 300 g for 1 min to pellet down worms, then the supernatants were removed. Worms were washed four more times with 1 mL M9 buffer. One hundred milliliters of 10 mM levamisole was added in M9 buffer to paralyze the worms for 5 min, followed by sterilization using 200 mL 4% bleach solution for 2 min. Worms were then washed two more times with M9 buffer to remove bleach and levamisole solutions. An aliquot of liquid volume was taken as a negative control to assess background residual live bacteria after the last wash. Three hundred microliters of the remaining worms was then lysed by adding 1.0-mm sterilized zirconia beads in a bullet blender homogenizer (Next Advance, USA) at full speed for 10 min to release live bacteria into the buffer.

Gut microbiome composition

Microbial community genomic DNA was extracted from worm lysate using TGuide S96 Magnetic Soil/Stool DNA Kit (Tiangen biochemical technology, China) according to the manufacturer’s instructions. The V3–V4 hypervariable region of the bacterial 16S rRNA gene was amplified with primer pairs 338F (5′-ACTCCTACGGGAGGCAGCA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′). PCR started with an initial denaturation step at 95 °C for 5 min; followed by 25 cycles of 95 °C for 30 s, 50 °C for 30 s, and 72 °C for 40 s; and one final elongation step at 72 °C for 7 min. Purified amplicons were sequenced on an Illumina NovaSeq 6000 platform (Illumina, USA) with 2 × 250 bp paired-end reads. Raw reads were quality-trimmed using Trimmomatic v0.33 [61]. Microbial community analysis of the 16S rRNA gene amplicon data was carried out using the QIIME2 environment v2022.2 with several available plugins [54]. Briefly, Demultiplexed paired-end reads were trimmed to remove primer sequences with cutadapt v2.7 [52], and read pairs were joined using VSEARCH [55]. Error correction, trimming (to a length of 430 nucleotides), and ASV clustering were performed using Deblur [89]. ASVs obtained by Deblur were used for a taxonomic assignment using the QIIME feature-classifier plugin against the manual training database based on V3–V4 region of 16S rRNA sequences of the HQbiome isolates. Further statistical analyses and graphical visualization of the ASVs data were performed in R. The DESeq2 R package was used to identify genera differentially abundant across different time points.

Availability of data and materials

All data generated or analyzed during this study are included in this published article, its supplementary information files and publicly available repositories. Scripts and data used for custom computational methods are available at Github (https://github.com/Erique29/Litoditis-marina-microbiome) and has been archived to Zenodo (https://zenodo.org/records/13762618). All raw genome data have been deposited in NCBI under BioProject accession number PRJNA922334. Other code and data of HQbiome isolates used in this study have been made available through the Figshare (https://doi.org/10.6084/m9.figshare.27024721.v4). All sequencing data of this study are available in Figshare: genome assemblies and annotations (https://doi.org/10.6084/m9.figshare.27024721.v4), 16S sequence data of environmental microbes (https://doi.org/10.6084/m9.figshare.27038002), and 16S sequence data of gut microbiota (https://doi.org/10.6084/m9.figshare.27038038). The genome of W. oceani could be available from NCBI (BioProject PRJNA326567 [90, 91]). Other data files analyzed during this study are included in this published article and its supporting information files. Supporting data values for n < 6 individual data values reported in the figures are detailed in the Additional file 6.

References

Gilbert SF, Sapp J, Tauber AI. A symbiotic view of life: we have never been individuals. Q Rev Biol. 2012;87(4):325–41.

Bruijning M, Henry LP, Forsberg SKG, Metcalf CJE, Ayroles JF. Natural selection for imprecise vertical transmission in host–microbiota systems. Nat Ecol Evol. 2022;6(1):77–87.

Wippel K, Tao K, Niu Y, Zgadzaj R, Kiel N, Guan R, et al. Host preference and invasiveness of commensal bacteria in the Lotus and Arabidopsis root microbiota. Nat Microbiol. 2021;6(9):1150–62.

Borrel G, Brugere JF, Gribaldo S, Schmitz RA, Moissl-Eichinger C. The host-associated archaeome. Nat Rev Microbiol. 2020;18(11):622–36.

Muller DB, Vogel C, Bai Y, Vorholt JA. The plant microbiota: systems-level insights and perspectives. Annu Rev Genet. 2016;50(1):211–34.

Vandenkoornhuyse P, Quaiser A, Duhamel M, Le Van A, Dufresne A. The importance of the microbiome of the plant holobiont. New Phytol. 2015;206(4):1196–206.

Baldassarre L, Reitzel AM, Fraune S. Genotype–environment interactions determine microbiota plasticity in the sea anemone Nematostella vectensis. PLos Biol. 2023;21(1):e3001726.

Phelps D, Brinkman NE, Keely SP, Anneken EM, Catron TR, Betancourt D, et al. Microbial colonization is required for normal neurobehavioral development in zebrafish. Sci Rep. 2017;7(1):11244.

Davis DJ, Bryda EC, Gillespie CH, Ericsson AC. Microbial modulation of behavior and stress responses in zebrafish larvae. Behav Brain Res. 2016;311:219–27.

Jaworska K, Koper M, Ufnal M. Gut microbiota and renin-angiotensin system: a complex interplay at local and systemic levels. Am J Physiol Gastrointest Liver Physiol. 2021;321(4):G355–66.

Fung TC, Olson CA, Hsiao EY. Interactions between the microbiota, immune and nervous systems in health and disease. Nat Neurosci. 2017;20(2):145–55.

de Vos WM, Tilg H, Van Hul M, Cani PD. Gut microbiome and health: mechanistic insights. Gut. 2022;71(5):1020–32.

Zhang L, Gualberto DG, Guo X, Correa P, Jee C, Garcia LR. TMC-1 attenuates C. elegans development and sexual behaviour in a chemically defined food environment. Nat Commun. 2015;6(1):6345.

Zhang F, Weckhorst JL, Assie A, Hosea C, Ayoub CA, Khodakova AS, et al. Natural genetic variation drives microbiome selection in the Caenorhabditis elegans gut. Curr Biol. 2021;31(12):2603–18.

Samuel BS, Rowedder H, Braendle C, Felix MA, Ruvkun G. Caenorhabditis elegans responses to bacteria from its natural habitats. Proc Natl Acad Sci U S A. 2016;113(27):E3941–9.

Qi B, Kniazeva M, Han M. A vitamin-B2-sensing mechanism that regulates gut protease activity to impact animal’s food behavior and growth. Elife. 2017;6: e26243.

Wei W, Ruvkun G. Lysosomal activity regulates Caenorhabditis elegans mitochondrial dynamics through vitamin B12 metabolism. Proc Natl Acad Sci U S A. 2020;117(33):19970–81.

Watson E, MacNeil LT, Ritter AD, Yilmaz LS, Rosebrock AP, Caudy AA, et al. Interspecies systems biology uncovers metabolites affecting C. elegans gene expression and life history traits. Cell. 2014;156(4):759–70.

Qi B, Han M. Microbial siderophore enterobactin promotes mitochondrial iron uptake and development of the host via interaction with ATP synthase. Cell. 2018;175(2):571–82.

Zhang J, Li X, Olmedo M, Holdorf AD, Shang Y, Artal-Sanz M, et al. A delicate balance between bacterial iron and reactive oxygen species supports optimal C. elegans development. Cell Host Microbe. 2019;26(3):400–11.

O’Donnell MP, Fox BW, Chao PH, Schroeder FC, Sengupta P. A neurotransmitter produced by gut bacteria modulates host sensory behaviour. Nature. 2020;583(7816):415–20.

Zimmermann J, Obeng N, Yang W, Pees B, Petersen C, Waschina S, et al. The functional repertoire contained within the native microbiota of the model nematode Caenorhabditis elegans. ISME J. 2020;14(1):26–38.

Berg M, Monnin D, Cho J, Nelson L, Crits-Christoph A, Shapira M. TGFbeta/BMP immune signaling affects abundance and function of C. elegans gut commensals. Nat Commun. 2019;10(1):604.

Dirksen P, Marsh SA, Braker I, Heitland N, Wagner S, Nakad R, et al. The native microbiome of the nematode Caenorhabditis elegans: gateway to a new host–microbiome model. BMC Biol. 2016;14(1):38.

Derycke S, De Meester N, Rigaux A, Creer S, Bik H, Thomas WK, et al. Coexisting cryptic species of the Litoditis marina complex (Nematoda) show differential resource use and have distinct microbiomes with high intraspecific variability. Mol Ecol. 2016;25(9):2093–110.

Xie Y, Zhang P, Xue B, Cao X, Ren X, Wang L, et al. Establishment of a marine nematode model for animal functional genomics, environmental adaptation and developmental evolution. BioRxiv. 2020. https://doi.org/10.1101/2020.03.06.980219.

Zhao L, Gao F, Gao S, Liang Y, Long H, Lv Z, et al. Biodiversity-based development and evolution: the emerging research systems in model and non-model organisms. Sci China Life Sci. 2021;64(8):1236–80.

Jungblut AD, Lovejoy C, Vincent WF. Global distribution of cyanobacterial ecotypes in the cold biosphere. ISME J. 2010;4(2):191–202.

Jung P, Briegel-Williams L, Schermer M, Budel B. Strong in combination: polyphasic approach enhances arguments for cold-assigned cyanobacterial endemism. Microbiologyopen. 2019;8(5): e00729.

Eilers H, Pernthaler J, Glockner FO, Amann R. Culturability and In situ abundance of pelagic bacteria from the North Sea. Appl Environ Microbiol. 2000;66(7):3044–51.

Galperin MY, Makarova KS, Wolf YI, Koonin EV. Expanded microbial genome coverage and improved protein family annotation in the COG database. Nucleic Acids Res. 2015;43(Database issue):D261–269.

Lieven C, Beber ME, Olivier BG, Bergmann FT, Ataman M, Babaei P, et al. MEMOTE for standardized genome-scale metabolic model testing. Nat Biotechnol. 2020;38(3):272–6.

Zarecki R, Oberhardt MA, Reshef L, Gophna U, Ruppin E. A novel nutritional predictor links microbial fastidiousness with lowered ubiquity, growth rate, and cooperativeness. PLoS Comp Biol. 2014;10(7): e1003726.

Sato K, Yoshiga T, Hasegawa K. Involvement of vitamin B6 biosynthesis pathways in the insecticidal activity of Photorhabdus luminescens. Appl Environ Microbiol. 2016;82(12):3546–53.

Datsenko KA, Wanner BL. One-step inactivation of chromosomal genes in Escherichia coli K-12 using PCR products. Proc Natl Acad Sci U S A. 2000;97(12):6640–5.

Miller DL, Roth MB. Hydrogen sulfide increases thermotolerance and lifespan in Caenorhabditis elegans. Proc Natl Acad Sci U S A. 2007;104(51):20618–22.

Joint I, Mühling M, Querellou J. Culturing marine bacteria–an essential prerequisite for biodiscovery. Microb Biotechnol. 2010;3(5):564–75.

Wang Y, Hekimi S. Understanding ubiquinone. Trends Cell Biol. 2016;26(5):367–78.

Nordman T, Xia L, Bjorkhem-Bergman L, Damdimopoulos A, Nalvarte I, Arner ES, et al. Regeneration of the antioxidant ubiquinol by lipoamide dehydrogenase, thioredoxin reductase and glutathione reductase. BioFactors. 2003;18(1–4):45–50.

Yang YY, Gangoiti JA, Sedensky MM, Morgan PG. The effect of different ubiquinones on lifespan in Caenorhabditis elegans. Mech Ageing Dev. 2009;130(6):370–6.

Hihi AK, Gao Y, Hekimi S. Ubiquinone is necessary for Caenorhabditis elegans development at mitochondrial and non-mitochondrial sites. J Biol Chem. 2002;277(3):2202–6.

Hihi AK, Kebir H, Hekimi S. Sensitivity of Caenorhabditis elegans clk-1 mutants to ubiquinone side-chain length reveals multiple ubiquinone-dependent processes. J Biol Chem. 2003;278(42):41013–8.

Jonassen T, Larsen PL, Clarke CF. A dietary source of coenzyme Q is essential for growth of long-lived Caenorhabditis elegans clk-1 mutants. Proc Natl Acad Sci U S A. 2001;98(2):421–6.

Jonassen T, Marbois BN, Faull KF, Clarke CF, Larsen PL. Development and fertility in Caenorhabditis elegans clk-1 mutants depend upon transport of dietary coenzyme Q8 to mitochondria. J Biol Chem. 2002;277(47):45020–7.

Severance S, Rajagopal A, Rao AU, Cerqueira GC, Mitreva M, El-Sayed NM, et al. Genome-wide analysis reveals novel genes essential for heme homeostasis in Caenorhabditis elegans. PLoS Genet. 2010;6(7): e1001044.

Read AD, Bentley RE, Archer SL, Dunham-Snary KJ. Mitochondrial iron–sulfur clusters: Structure, function, and an emerging role in vascular biology. Redox Biol. 2021;47: 102164.

Correia-Melo C, Kamrad S, Tengölics R, Messner CB, Trebulle P, Townsend S, et al. Cell-cell metabolite exchange creates a pro-survival metabolic environment that extends lifespan. Cell. 2023;186(1):63–79.

Snell TW, Johnston RK. Glycerol extends lifespan of Brachionus manjavacas (Rotifera) and protects against stressors. Exp Gerontol. 2014;57:47–56.

Lee S-J, Murphy CT, Kenyon C. Glucose shortens the life span of C. elegans by downregulating DAF-16/FOXO activity and aquaporin gene expression. Cell Metab. 2009;10(5):379–91.

Ghaddar A, Mony VK, Mishra S, Berhanu S, Johnson JC, Enriquez-Hesles E, et al. Increased alcohol dehydrogenase 1 activity promotes longevity. Curr Biol. 2023;33(6):1036–46.

Holterman M, Schratzberger M, Helder J. Nematodes as evolutionary commuters between marine, freshwater and terrestrial habitats. Biol J Linn Soc. 2019;128(3):756–67.

Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet j. 2011;17(1):10–2.

Callahan BJ, Wong J, Heiner C, Oh S, Theriot CM, Gulati AS, et al. High-throughput amplicon sequencing of the full-length 16S rRNA gene with single-nucleotide resolution. Nucleic Acids Res. 2019;47(18): e103.

Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019;37(8):852–7.

Rognes T, Flouri T, Nichols B, Quince C, Mahe F. VSEARCH: a versatile open source tool for metagenomics. PeerJ. 2016;4: e2584.

Douglas GM, Maffei VJ, Zaneveld JR, Yurgel SN, Brown JR, Taylor CM, et al. PICRUSt2 for prediction of metagenome functions. Nat Biotechnol. 2020;38(6):685–8.

Caspi R, Billington R, Keseler IM, Kothari A, Krummenacker M, Midford PE, et al. The MetaCyc database of metabolic pathways and enzymes–a 2019 update. Nucleic Acids Res. 2020;48(D1):D445–53.

Hartigan JA, Wong MA. Algorithm AS 136: A k-means clustering algorithm. J Roy Stat Soc Ser C (Appl Stat). 1979;28(1):100–8.

Xie Y, Zhang L. Transcriptomic and proteomic analysis of marine nematode Litoditis marina acclimated to different salinities. Genes (Basel). 2022;13(4): 651.

Lim HJ, Lee EH, Yoon Y, Chua B, Son A. Portable lysis apparatus for rapid single-step DNA extraction of Bacillus subtilis. J Appl Microbiol. 2016;120(2):379–87.

Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biol. 2012;19(5):455–77.

Li H. Minimap and miniasm: fast mapping and de novo assembly for noisy long sequences. Bioinformatics. 2016;32(14):2103–10.

Vaser R, Sovic I, Nagarajan N, Sikic M. Fast and accurate de novo genome assembly from long uncorrected reads. Genome Res. 2017;27(5):737–46.

Wick RR, Judd LM, Gorrie CL, Holt KE. Unicycler: Resolving bacterial genome assemblies from short and long sequencing reads. PLoS Comp Biol. 2017;13(6): e1005595.

Chin CS, Alexander DH, Marks P, Klammer AA, Drake J, Heiner C, et al. Nonhybrid, finished microbial genome assemblies from long-read SMRT sequencing data. Nat Methods. 2013;10(6):563–9.

Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60.

Simao FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31(19):3210–2.

Gurevich A, Saveliev V, Vyahhi N, Tesler G. QUAST: quality assessment tool for genome assemblies. Bioinformatics. 2013;29(8):1072–5.

Alikhan NF, Petty NK, Ben Zakour NL, Beatson SA. BLAST ring image generator (BRIG): simple prokaryote genome comparisons. BMC Genomics. 2011;12(1): 402.

Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics. 2014;30(14):2068–9.

Huerta-Cepas J, Forslund K, Coelho LP, Szklarczyk D, Jensen LJ, von Mering C, et al. Fast genome-wide functional annotation through orthology assignment by eggNOG-mapper. Mol Biol Evol. 2017;34(8):2115–22.

Jain C, Rodriguez RL, Phillippy AM, Konstantinidis KT, Aluru S. High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat Commun. 2018;9(1):5114.

Wu Y-W. ezTree: an automated pipeline for identifying phylogenetic marker genes and inferring evolutionary relationships among uncultivated prokaryotic draft genomes. BMC Genomics. 2018;19(1):7–16.

Letunic I, Bork P. Interactive Tree Of Life (iTOL) v5: an online tool for phylogenetic tree display and annotation. Nucleic Acids Res. 2021;49(W1):W293–6.

Revell LJ. phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol Evol. 2012;3(2):217–23.

Jombart T, Balloux F, Dray S. adephylo: new tools for investigating the phylogenetic signal in biological traits. Bioinformatics. 2010;26(15):1907–9.

Zimmermann J, Kaleta C, Waschina S. gapseq: informed prediction of bacterial metabolic pathways and reconstruction of accurate metabolic models. Genome Biol. 2021;22(1):81.

Saier MH Jr, Reddy VS, Tamang DG, Vastermark A. The transporter classification database. Nucleic Acids Res. 2014;42(database issue):d251–258.

Holzhutter HG. The principle of flux minimization and its application to estimate stationary fluxes in metabolic networks. Eur J Biochem. 2004;271(14):2905–22.

Dirksen P, Assié A, Zimmermann J, Zhang F, Tietje A-M, Marsh SA, et al. CeMbio–the caenorhabditis elegans microbiome resource. G3. 2020;10(9):3025–39.

Mahadevan R, Schilling CH. The effects of alternate optimal solutions in constraint-based genome-scale metabolic models. Metab Eng. 2003;5(4):264–76.

Bauer E, Zimmermann J, Baldini F, Thiele I, Kaleta C. BacArena: Individual-based metabolic modeling of heterogeneous microbes in complex communities. PLoS Comp Biol. 2017;13(5): e1005544.

Lagesen K, Hallin P, Rodland EA, Staerfeldt HH, Rognes T, Ussery DW. RNAmmer: consistent and rapid annotation of ribosomal RNA genes. Nucleic Acids Res. 2007;35(9):3100–8.

Suzuki R, Shimodaira H. Pvclust: an R package for assessing the uncertainty in hierarchical clustering. Bioinformatics. 2006;22(12):1540–2.

Arkin AP, Cottingham RW, Henry CS, Harris NL, Stevens RL, Maslov S, et al. KBase: the United States department of energy systems biology knowledgebase. Nat Biotechnol. 2018;36(7):566–9.

Cáceres MD, Legendre P. Associations between species and groups of sites: indices and statistical inference. Ecology. 2009;90(12):3566–74.

Genuer R, Poggi JM, Tuleau-Malot C. VSURF: an R package for variable selection using random forests. R J. 2015;7(2):19–33.

Amir A, McDonald D, Navas-Molina JA, Kopylova E, Morton JT, Zech XuZ, et al. Deblur rapidly resolves single-nucleotide community sequence patterns. mSystems. 2017;2(2):e00191–00116.

Du ZJ, Wang ZJ, Zhao JX, Chen GJ. Woeseia oceani gen. nov., sp. nov., a chemoheterotrophic member of the order Chromatiales, and proposal of Woeseiaceae fam. nov. Int J Syst Evol Microbiol. 2016;66(1):107–12.

Du ZJ, Wang ZJ, Zhao JX, Chen GJ. Woeseia oceani Genome sequencing. GenBank https://www.ncbi.nlm.nih.gov/nuccore/1041525582 (2016).