Crab Survival and Growth
There was no mortality in dietary treatments of FO and LOCO throughout the feeding trial, while 77% survival occurred in diet LOCOAB treatment (Table S3). Diet LOCOAB also resulted in lower percentage of molting and specific growth rate as compared to diets FO and LOCO. This also resulted in higher feed conversion ratio value in LOCOAB treatment, despite similar feed intake between crabs from different dietary treatments.
Overview Of The Metagenomic Sequencing Data
A total of 559 million paired-end reads, with an average of 47 million reads per sample was generated from shotgun metagenomic sequencing of the WT, FO, LOCO, and LOCOAB samples (Table 1). After trimming, 537 million clean reads were assembled into 52,275-85,353 contigs with a total length of 94-174 Mb and N50 contig length of 2,615-5,843 bp. A non-redundant protein-coding gene set containing ~1.89 million ORFs was predicted from the 12 metagenomes.
Table 1
Summary of the shotgun metagenomic sequencing data of S. olivacea gut microbiota from wild-caught crabs (WT), crabs fed diet FO, LOCO, or LOCOAB.
Samples | No. of raw reads | No. of clean reads | No. of contigs | N50 (bp) | Assembly size (Mb) | No. of ORFs | Average ORF length (bp) |
WT1 | 43,931,876 | 42,488,724 | 53,513 | 2,615 | 95.30 | 126,021 | 665 |
WT2 | 48,522,778 | 46,858,652 | 56,482 | 3,018 | 105.76 | 137,920 | 673 |
WT3 | 49,736,880 | 48,278,740 | 52,275 | 2,815 | 93.95 | 124,356 | 666 |
FO1 | 54,927,728 | 52,825,174 | 85,353 | 5,395 | 173.66 | 218,756 | 699 |
FO2 | 44,267,802 | 42,378,718 | 71,426 | 5,843 | 156.65 | 193,531 | 716 |
FO3 | 40,904,536 | 39,423,340 | 64,900 | 5,791 | 147.33 | 180,376 | 723 |
LOCO1 | 48,501,854 | 46,462,276 | 66,117 | 4,796 | 133.61 | 169,367 | 702 |
LOCO2 | 45,912,916 | 44,165,248 | 61,136 | 4,742 | 126.62 | 159,463 | 709 |
LOCO3 | 47,387,810 | 45,414,300 | 65,324 | 4,939 | 132.59 | 167,666 | 705 |
LOCOAB1 | 49,051,030 | 46,775,560 | 63,356 | 2,894 | 110.96 | 142,647 | 676 |
LOCOAB2 | 42,311,286 | 40,324,828 | 60,475 | 2,642 | 102.10 | 133,164 | 669 |
LOCOAB3 | 43,761,230 | 41,718,530 | 58,772 | 2,824 | 101.76 | 131,813 | 672 |
The 16S rRNA amplicon sequencing produced 2.95 million reads from the 12 individual samples, ranging from 212,812 to 292,709 reads for each sample (Table 2). These sequences were delineated into 18,514 OTUs, corresponding to an average of 1,543 OTUs per sample. For community richness and diversity comparison, alpha diversity parameters were calculated from the proportion of OTUs. The Chao 1 index, an indicator of microbiota community richness, varied from 1,410 to 2,124. Among the different treatments, FO and LOCO samples exhibited higher OTUs and Chao1 index values than WT and LOCOAB. Similarly, Shannon index value was also highest in FO. The Good's coverage estimator of sequencing completeness ranged from 0.998 to 0.999, indicating high species coverage within the samples. Additionally, the species accumulation curve appeared to flatten after 8, indicating that the number of samples was sufficient to reflect the species abundance (Fig. S2A). LOCO samples showed the widest and smoothest curve in rank abundance analysis, indicating the highest species richness and, concomitantly, the species uniformity among sample groups (Fig. S2B). Species richness and community evenness were lowest in the LOCOAB samples.
Table 2
Summary of the 16S rRNA sequencing data and alpha-diversity indexes of S. olivacea gut microbiota from wild-caught crabs (WT), crabs fed diet FO, LOCO, or LOCOAB.
Samples | No. of reads | OTU | Shannon | Simpson | Chao1 | Good's coverage |
WT1 | 292,709 | 1,489 | 2.85 | 0.24 | 1,820 | 0.999 |
WT2 | 268,093 | 1,482 | 3.03 | 0.19 | 1,822 | 0.999 |
WT3 | 220,682 | 1,332 | 2.75 | 0.25 | 1,709 | 0.998 |
FO1 | 212,812 | 1,648 | 3.92 | 0.05 | 1,977 | 0.998 |
FO2 | 244,073 | 1,746 | 3.90 | 0.06 | 2,019 | 0.998 |
FO3 | 262,948 | 1,815 | 4.04 | 0.05 | 2,107 | 0.999 |
LOCO1 | 214,435 | 1,752 | 3.80 | 0.08 | 2,053 | 0.998 |
LOCO2 | 241,125 | 1,800 | 3.45 | 0.14 | 2,124 | 0.998 |
LOCO3 | 247,755 | 1,874 | 3.75 | 0.10 | 2,123 | 0.999 |
LOCOAB1 | 264,262 | 1,214 | 3.38 | 0.10 | 1,440 | 0.999 |
LOCOAB2 | 217,363 | 1,144 | 3.29 | 0.11 | 1,410 | 0.999 |
LOCOAB3 | 261,232 | 1,218 | 3.42 | 0.09 | 1,462 | 0.999 |
The PCoA and NMDS analyses of beta-diversity were performed to visualize the differences in bacterial communities among samples (Fig. 1A and 1B). Principal component 1 (PC1) and PC2 accounted for 72.34% and 17.38% of the composition variance, respectively, reflecting the dissimilarity in the bacterial community composition among different dietary groups. Additionally, the hierarchical clustering tree showed a clear distinction among these groups while the replicates within each group are consistently grouped (Fig. 1C). The Venn diagram showed 684 shared OTUs among the four groups and 389, 133, 165, and 153 unique OTUs within the WT, FO, LOCO, and LOCOAB groups, respectively (Fig. 1D). Treatments FO and LOCO shared the highest number of OTUs (1,762), followed by WT and FO (1,245), WT and LOCO (1,231), and lastly, LOCO and LOCOAB (1,166).
Gut Microbiota Taxonomic Composition
At the phylum level, Proteobacteria (21.8-57.4%), Firmicutes (11.8-38.8%), Bacteroidota (2.7-32.0%), Fusobacteriota (1.3-11.1%), and Tenericutes (0.3-52.2%) were the core groups in all types of samples, accounting for 86.2% of the total reads (Fig. 2A and 2B). The phylum Tenericutes was reclassified into the Bacilli class of Firmicutes following the SILVA database and therefore, not annotated in the 16S rRNA amplicon analysis. The WT microbiota samples were dominated by Firmicutes and Tenericutes, while in FO samples, Proteobacteria and Fusobacteria were the most abundant. Proteobacteria also dominated the LOCO microbiota, while for diet LOCOAB treatment, Bacteroidota, Proteobacteria, and Spirochaetota were most abundant. Experimental diets treatment reduced the abundance of Firmicutes and Tenericutes, increasing the proportion of Proteobacteria, Fusobacteria, Bacteroidota, Spirochaetota, and Campylobacterota, respectively. The microbiota of LOCO treatment showed a higher Proteobacteria, Firmicutes, and Tenericutes, and lower Bacteroidota, Fusobacteria, Spirochaetota, and Campylobacterota abundance as compared to diet FO. Treatment with diet LOCOAB lowered the abundance of Proteobacteria, Firmicutes, Tenericutes, and Fusobacteria, while Bacteroidota and Spirochaetota increased.
The Proteobacteria phylum was mainly represented by Gammaproteobacteria, predominantly from the Vibrionaceae family (Fig. S3 and S4). Fusobacteria in the FO microbiota was represented by the Fusobacteriaceae family from the Fusobacteriia class. In WT samples, the Entomoplasmatales Incertae Sedis family of Bacilli class was predominant in Firmicutes, and the Mycoplasmataceae family of Mollicutes class was most abundant among Tenericutes. Within the Bacteroidota, the Marinilabilaceae and Prolixibacteraceae families dominated the LOCOAB microbiota. At the genus level, both the 16S and shotgun datasets showed the dominance of Candidatus Hepatoplasma in the WT microbiota (Fig. 2C and 2D). Compared to animals fed experimental diets, there was a reduction in Candidatus Hepatoplasma, alongside Mycoplasma, Photobacterium, Spiroplasma, Paraclostridium, Bacillus, and Hypnocyclicus in the microbiota of experimentally fed animals. In contrast, the abundance of Vibrio, Sediminispirochaeta, Sunxiuqinia, Ruegeria, Arcobacter, Malaciobacter, and Carboxylicivirga was elevated. In diet LOCO, Vibrio, Candidatus Hepatoplasma, Shewanella, and Ferrimonas showed higher abundance than diet FO, while the abundance of Sediminispirochaeta, Carboxylicivirga, Sunxiuginia, Psychrilyobacter, Arcobacter, Malaciobacter, Propionigenium, and Halarcobacter was lower. For a more in-depth comparison between diet LOCO and LOCOAB, the Statistical Analysis of Metagenome Profiles analysis was carried out (Fig. S5). In LOCOAB samples, genera such as Vibrio, Candidatus Hepatoplasma, Ruegeria, Shewanella Mycoplasma, Spiroplasma, Photobacterium, and Psychrilyobacter were reduced, while Sediminispirochaeta, Carboxylicivirga, Sunxiuqinia, and Oceanispirochaeta were elevated.
Supervised comparisons by LEfSe (LDA > 4.0) were also conducted to provide an overview of differences in taxonomic profiles for all dietary treatments (Fig. 3). The analysis identified 50 differentially abundant taxa ranging from the phylum to OTU level (12 in WT; 16 in FO; 7 in LOCO; 15 in LOCOAB).
S. olivacea Gut Microbiota Functions
A total of 84.4% of the unique gene sets from all microbiota samples were assigned into functional groups (Table S4). These genes were further annotated to 44 pathways at KEGG level 2, with the highest enrichment in global and overview maps (25.4%), followed by carbohydrate (8.3%) and energy (5.2%) metabolisms (Fig. S6). The relative composition of these three categories at KEGG level 3 is illustrated in Fig. S7.
A heatmap with normalized values of functional abundance was illustrated to compare the enriched level 2 categories between the four sample groups (Fig. 4A). The WT crab microbiota showed a higher relative abundance of nucleotide metabolism, transcription, translation, folding, sorting and degradation, replication and repair, and signaling molecules and interaction pathways. The microbiota of crabs fed experimental diets were enriched in global and overview maps, numerous metabolisms including lipid, amino acid, other amino acids, cofactors and vitamins, terpenoids and polyketides, and biosynthesis of secondary metabolites. Among the experimental diets, LOCO resulted in enrichment in cell motility, cellular community-prokaryotes, membrane transport, signal transduction, and xenobiotics degradation and metabolism. For diet LOCOAB, highest abundance of global and overview maps, and metabolisms of major nutrients were obtained. In addition, the metagenome gene sets were also mapped against the Comprehensive Antibiotic Resistance Database (CARD) to determine the occurrence of antibiotic resistance genes (ARGs) in the S. olivacea gut microbiota (Fig. S8). We recovered ARGs from 15 main families, corresponding to seven drug mechanisms in CARD.
In relevance to the main objective of this study, a heatmap illustration to compare the KEGG Level 3 categories within lipid metabolism between dietary treatments was built (Fig. 4B). Feeding crabs with experimental diets enriched various lipid metabolism pathways, including fatty acid biosynthesis and fatty acid degradation. Compared to diet FO, pathways on fatty acid degradation, synthesis and degradation of ketone bodies, glycerophospholipid metabolism, ether lipid metabolism, linoleic acid metabolism, arachidonic acid metabolism, and biosynthesis of unsaturated fatty acids were enriched in the LOCO microbiota. Concomitantly, diet LOCOAB lowered the abundance of these functions while causing an increase in sphingolipid metabolism, steroid hormone biosynthesis, fatty acid elongation, and fatty acid biosynthesis.
Within lipid metabolism, the fatty acid biosynthesis and biosynthesis of unsaturated fatty acids are key pathways towards the biosynthesis of LC-PUFA. Numerous genes encoding for the enzymes are present in both pathways (Fig. S9 and S10). There was an increase in the abundance of several sequences assigned to the KEGG function biosynthesis of unsaturated fatty acid in diet LOCO treatment. Among these are acyl-CoA oxidase, stearoyl-CoA desaturase, and acyl-coenzyme A thioesterases, which were also impeded in diet LOCOAB (Table S5). Conversely, there was higher abundance of fatty acid biosynthesis genes in the microbiota of LOCOAB treatment. The complete pathways for acetate and butyrate formation and a partial propionate pathway are present in S. olivacea gut microbiota (Fig. S11-S13). Heatmap analysis showed the LOCOAB treatment increased the function of short-chain fatty acids (SCFAs) production as compared to LOCO treatment (Fig. 4C).
In terms of bacteria taxa, genes involved in lipid metabolism were mostly contributed by Proteobacteria and represented by the families Vibrionaceae, Rhodobacteraceae, and Shewanellacea (Fig. 5A). At the genus level, lipid metabolism in WT microbiota was contributed by Vibrio, Shewanella, and Clostridium (Fig. 5B). In diet LOCO, the Ruegeria and Shewanella genera contributed higher towards lipid metabolism as compared to diet FO. While diet LOCOAB did not affect the contribution levels of these two genera, there was an increase in contribution from Vibrio, Arcobacter, Clostridium, and Bacteroides.
Several marine species, in particular Shewanella and Vibrio, can synthesize PUFA de novo through an anaerobic route, using the fatty acid synthase/polyketide synthase (FAS/PKS)-like enzyme system [40, 41]. The gene cluster encoding this system consists of four ORFs, represented by the pfaABCD genes [42]. Since this pathway is not included in the KEGG database, we examined the S. olivacea gut shotgun metagenomic data sets for the diversity and distribution of the keto-acyl synthase (KS) domain harbored within the pfaA homolog. The KS sequences recovered in this study form a large monophyletic clade with known KSs from Vibrio, Shewanella, and Photobacterium, showing conservation of the KS domain across species (Fig. 6). A total of 362 KS sequences were recovered, with 36, 135, 127, and 64 sequences from WT, FO, LOCO, and LOCOAB samples, respectively. Therefore, in contrast to SCFAs pathways, the inclusion of oxolinic acid impeded the PKS pathway in S. olivacea.
In hepatopancreas, the highest levels of DHA were obtained with diet FO (Table S6). Diet LOCO resulted in the deposition of α-linolenic acid, ALA and linoleic acid, LA, and also some known intermediates of the LC-PUFA biosynthesis pathway such as 18:4n-3 and 22:5n-3, which implies some degree of LC-PUFA biosynthesis. Interestingly, despite having the same lipid level and fatty acids profile, the percentage of ALA and LA were significantly higher in diet LOCOAB. Crabs fed diet FO also have the highest deposition of all three LC-PUFA in muscle tissue. There was no significant difference in LC-PUFA levels between diet LOCO and LOCOAB.