Experimental design. Fish were bred and raised as described in the section ‘Zebrafish and diets’. At 40 dpf (juvenile stage) they were fed diet A for 1 week for acclimatisation to dry feed pellets. At 47 dpf, fish were randomly distributed into tanks and fed one of the diets (A, B or C) for 3 weeks. Survival and growth were measured before and during the whole experiment. After 1 week (54 dpf) and after 3 weeks (68 dpf) of feeding the fish. Gut samples were collected for histological, metatranscriptomic and microbiome analyses

Combinatorial approach employed: total RNA was extracted from single zebrafish gut fed on different diets for both timepoints. Aliquots of total RNA were used for cDNA. For the 16S rRNA gene profiling, amplicon libraries of the V4 region of the 16S RNA gene were generated from the cDNA synthetized. NG-Tax 2.0 Galaxy was sued to obtain the ASVs. Several packages of R v4.1.2., Canoco v5.15 and Cytoscape v3.9.1 were used for results visualization. For transcriptomics, the cDNA libraries were sent to NovaSeq 6000 PE150 for sequencing. MetaPhlAn 3.0 (Beghini et al. 2021) and KneadData were used to trim the overrepresented sequences. Nf-core/rnaseq Nextflow pipeline was used for processing of the reads with the GRCz11 genome assembly. The results were visualized by R v.4.1.2, Canoco v5.15 and ErmineJ was used for the GO Enrichment analysis. The histological samples were extracted and embedded in paraffin and sectioned using a microtome. AB-PAS and HIC stains were automated. Samples were digitally scanned and an automated quantification of the histological parameters was performed using VIS v.2019.07 and Canoco v5.15 and GraphPad Prism v9.0.0 to visualize the results. The data integration was performed using heatmaps of normalized relevant parameters from all datasets, both timepoints and all diets

Alpha-diversity indexes for richness (observed ASVs and Chao1) and for diversity (Shannon, Inverse Simpson, Fisher and Phylogenetic Diversity). *p ≤ 0.05, Ordinary one-way ANOVA after confirming normally distributed data by Shapiro–Wilk test. Whiskers: min. to max. shall all points with median

A Principal component analysis exploring the interaction of diet and time. The x axis separates the samples after 1 week feeding (54 dpf) from samples after 3 weeks feeding (68 dpf) and explains 16.87% of the variation observed. B Redundancy analysis of samples after 3 weeks of feeding (68 dpf), the x axis separates saponin from butyrate fed fish and explains 7.94% of the microbial differences observed and the y axis separates the control from the saponin fed fish and explains 3.45% of the microbial differences observed. The microbial communities changed significantly due to diets (p = 0.018). The relative abundance of the most discriminative genera are depicted with boxplots around the RDA. In both analyses, the top 15 most distinctive genera are represented with black arrows. The direction of the arrows correlated with the dietary treatments and the timepoints and their length correlate with the strength of the correlation. **p ≤ 0.01, *** p ≤ 0.005 one-way ANOVA test or Kruskal–Wallis test after testing for normality on data distribution by Shapiro–Wilk test. No false discovery rate performed. Whiskers: min. to max. shall all points with median

Taxa connectivity: taxa included when prevalence is ≥ in 3/10 samples, abundance is ≥ 10 counts in 1 M and significance ≤ 0.1. The lines inform about the nature of the taxa interaction: the thickness of the lines represents the strength of the correlation (r-score value) and the shape of the lines represents the direction of the correlation, straight lines mean positive correlation (co-occurrence) whereas dashed lines mean negative correlation of the pairs of taxa (anti-occurrence). A Pairs of taxa co- and anti-occurring for all the diets: in black the interactions occurring in all three diets whereas in grey the interactions not occurring in all diets. Node size corresponds to average abundance of taxa for all diets at 68 dpf. B In red the interactions occurring in the control fed fish and not in the other two diets. Node size corresponds to average abundance of taxa for control diet at 68 dpf. C In green the interactions occurring in the butyrate fed fish and not in the other two diets. Node size corresponds to average abundance of taxa for butyrate diet at 68 dpf. D In yellow the interactions occurring in the saponin fed fish and not in the other two diets. Node size corresponds to average abundance of taxa for saponin diet at 68 dpf

Effects of butyrate and saponin on the host gut transcriptome. A Bray–Curtis distances to examine the dissimilarity of the host transcriptome across diets and timepoints. ** p ≤ 0.01, Kruskal–Wallis test after testing for non-normally distributed data by Shapiro–Wilk test. Whiskers: min. to max. shall all points with median. B Network depicting transcriptomic regulation of butyrate and saponin supplemented diets vs control diet at 68 dpf. Each node is a GO term and the node border represent the log2 fold-change of the control diet vs the saponin supplemented diet and the node fill represent the log2 fold-change of the control diet vs the butyrate supplemented diet. The edges connect nodes containing at least 10 genes and sharing 50% of the contained genes. Related GO terms are encircled encompassing canonical pathways. Shared effects on the gut transcriptome can be observed when edge and fill of a node have the same color in the network: up-regulation -in red- and down-regulation -in blue- compared to the control feed. C Immune response-associated GO terms and particularly inflammatory response analysed in fish fed a control, butyrate and saponin diet. Genes are expressed in tpm and scaled colored per individual gene value. The heatmap contained genes color-scaled per individual gene that reflect the individual within group fish-to-fish variation

Effects of butyrate and saponin on the host gut transcriptome. A Bray–Curtis distances to examine the dissimilarity of the host transcriptome across diets and timepoints. ** p ≤ 0.01, Kruskal–Wallis test after testing for non-normally distributed data by Shapiro–Wilk test. Whiskers: min. to max. shall all points with median. B Network depicting transcriptomic regulation of butyrate and saponin supplemented diets vs control diet at 68 dpf. Each node is a GO term and the node border represent the log2 fold-change of the control diet vs the saponin supplemented diet and the node fill represent the log2 fold-change of the control diet vs the butyrate supplemented diet. The edges connect nodes containing at least 10 genes and sharing 50% of the contained genes. Related GO terms are encircled encompassing canonical pathways. Shared effects on the gut transcriptome can be observed when edge and fill of a node have the same color in the network: up-regulation -in red- and down-regulation -in blue- compared to the control feed. C Immune response-associated GO terms and particularly inflammatory response analysed in fish fed a control, butyrate and saponin diet. Genes are expressed in tpm and scaled colored per individual gene value. The heatmap contained genes color-scaled per individual gene that reflect the individual within group fish-to-fish variation

High-throughput quantitative histological analysis. A Representative pictures of all cell-types analyzed for all the diets and timepoints with cells of interest in dashed black lines per each group. B Redundancy analysis to examine the effect of time on the histological parameters analyzed. The x axis separated the samples by timepoints and explained 5.76% of the variation observed. The link of time and variation of the histological parameters was not significant (p = 0.066) C Redundancy analysis to examine the effect of diet on the histological parameters analyzed. The x axis separated the samples of butyrate fed fish from saponin and control fed fish and explained 5.86% of the variation explained. The y axis separated saponin fed fish from control fed fish and explained 5.07% of the variation observed. The link of time and variation of the histological parameters was significant (p = 0.036). The top 10 most distinctive histological parameters are depicted in black arrows. The direction of the arrows correlate with the dietary intervention and the length of the arrows represents the strength of the correlation. The boxplots around the RDA depicted the absorptive capacity and the percentage area of cells of interest compared to the total gut area per diet and timepoint. *p ≤ 0.05, one-way ANOVA test or Kruskal–Wallis test after testing for normality on data distribution by Shapiro–Wilk test. No false discovery rate performed. Whiskers: min. to max. shall all points with median

High-throughput quantitative histological analysis. A Representative pictures of all cell-types analyzed for all the diets and timepoints with cells of interest in dashed black lines per each group. B Redundancy analysis to examine the effect of time on the histological parameters analyzed. The x axis separated the samples by timepoints and explained 5.76% of the variation observed. The link of time and variation of the histological parameters was not significant (p = 0.066) C Redundancy analysis to examine the effect of diet on the histological parameters analyzed. The x axis separated the samples of butyrate fed fish from saponin and control fed fish and explained 5.86% of the variation explained. The y axis separated saponin fed fish from control fed fish and explained 5.07% of the variation observed. The link of time and variation of the histological parameters was significant (p = 0.036). The top 10 most distinctive histological parameters are depicted in black arrows. The direction of the arrows correlate with the dietary intervention and the length of the arrows represents the strength of the correlation. The boxplots around the RDA depicted the absorptive capacity and the percentage area of cells of interest compared to the total gut area per diet and timepoint. *p ≤ 0.05, one-way ANOVA test or Kruskal–Wallis test after testing for normality on data distribution by Shapiro–Wilk test. No false discovery rate performed. Whiskers: min. to max. shall all points with median

The heatmap brings together the main observations of each analysis and compare them per diet and timepoint. The more representative genera are illustrated with the average relative abundance per timepoint and diet. The taxa connectivity contained the amount of pairs of taxa that correlate to each other in a significant fashion ( p ≤ 0.05). The GO terms contain the transcripts per million (tpm) of all genes expressed in the dataset that collapsed under that GO term. All histological parameters are normalized and scaled (from 0 to 1). Each individual feature within the heatmap is normalized and colored from red (more present) to white (absent)

A Representative pictures of the fluorescent in vivo imaging of the gut area of Tg(mpeg1:mCherry / mpx:eGFPi¹¹⁴) larvae exposed to control media, 0.005 mg/ml and 0.01 mg/ml butyrate and 0.5 mg/ml and 0.7 mg/ml saponin. B Quantification of neutrophils and macrophages in the gut area of the zebrafish larvae (n = 10 in all groups except 0.7 mg/ml saponin where n = 3). *p ≤ 0.5, **** p ≤ 0.0001 Kruskal–Wallis test after testing for non-normally distributed data by Shapiro–Wilk test. Whiskers: min. to max. shall all points with median