Serum was extracted from blood samples by centrifugation and stored at −80°C. Concentrations of liver function variables in serum, i.e., gamma glutamyl transferase (GGT), total bilirubin (TB), total bile acid (TBA), albumin (ALB), aspartate aminotransferase (AST), alanine aminotransferase (ALT), and alkaline phosphatase (ALP), were measured by an automatic biochemical analyzer (Roche Diagnostics, Ottweiler, Germany) according to the manufacturer’s instructions.

The concentrations of 23 cytokines in the serum samples were measured using a Bio-Plex Pro™ Rat Cytokine 23-Plex Assay kit (Bio-Rad Ltd., Hercules, CA, USA) as per the protocol of the manufacturer. These cytokines included macrophage inflammatory protein (MIP)-1α, MIP-3α, macrophage colony-stimulating factor (M-CSF), granulocyte colony-stimulating factor (G-CSF), granulocyte–macrophage colony-stimulating factor (GM-CSF), interferon gamma (IFN-γ), tumor necrosis factor alpha (TNF-α), interleukin (IL)-1α, IL-1β, IL-2, IL-4, IL-5, IL-6, IL-7, IL-10, IL-12p70, IL-13, IL-17A, IL-18, monocyte chemoattractant protein-1 (MCP-1), growth-related oncogene (GRO)/keratinocyte chemoattractant (KC), regulated upon activation, normal T cell expressed, and secreted (RANTES), and vascular endothelial growth factor (VEGF).

The liver tissue from the left liver lobe of each rat was dissected and fixed in 10% formalin solution, before being dehydrated and processed in paraffin using standard histological methods. The liver samples were stained and mounted on microscope slides. The liver injury severity was evaluated by a professional pathologist based on the Ishak scoring system (13).

DNA was extracted from the cecal samples by using a DNeasy PowerSoil kit (MoBio Laboratories Inc., Carlsbad, CA, USA). The extracted DNA was respectively amplified with bacterial primers (i.e., 341F/785R) and fungal primers (i.e., ITS3F/ITS4R) (14, 15). The PCR products were purified by using a DNA Clean and Concentrator Kit (Zymo Research, Irvine, CA, USA), and their concentrations were measured by using a Qubit™ dsDNA HS Assay Kit (Thermo Fisher Scientific Inc., Waltham, MA, USA). The purified PCR products were submitted for sequencing on Illumina NovaSeq 6000 platform (Illumina Inc. USA).

The sequencing data were processed with standard bioinformatic procedures, e.g., merge, quality filter, singleton removal, and chimera check. The sequences were clustered into groups of amplicon sequence variants (ASVs). Bacterial phylotypes were classified using QIIME2 against the Silva 138 database, while fungal phylotypes were classified against the UNITE fungal database. One rat in the LI10 cohort was not recruited in the current study, as it did not provide enough sequencing reads. All the other rats were recruited for the subsequent analyses.

One-way ANOVA tests and Mann–Whitney tests were used to determine the differences in 23 cytokines in LI09, LI10, PC, and NC cohorts. The cytokines with significant differences among the four cohorts were selected and transformed in log₂(X+1) for the cytokine profile analyses.

Partition around medoids (PAM) clustering analysis was performed to cluster the cytokine profiles of all the four cohorts, after an average silhouette analysis being conducted to determine the optimal number of clusters. The four cohorts were determined with two distinct cytokine profiles, i.e., CP1 and CP2. The rats in LI09 and LI10 cohorts with two different cytokine profiles were defined as CP1_LI09, CP2_LI09, CP1_LI10, and CP2_LI10 cohorts.

T-tests and Mann–Whitney tests were used to compare the liver function variables in CP1_LI09 and CP2_LI09 cohorts. The same tests were performed to compare the liver function variables in CP1_LI10 and CP2_LI10 cohorts.

Mann–Whitney tests were carried out to compare the Ishak scores of CP1_LI09 and CP2_LI09 cohorts and those of CP1_LI10 and CP2_LI10 cohorts.

The alpha diversity indices (i.e., observed species and Shannon and Pielou indices) of the bacterial and fungal microbiome in the LI09 cohorts with different cytokine profiles, and their control cohorts (i.e., PC and NC cohorts), were all calculated. One-way ANOVA was used to compare the alpha diversity indices of the four cohorts, and t-tests were performed for the comparisons of LI09 and control cohorts. The same strategy was carried out for the same microbiome composition comparisons of LI10 cohorts with different cytokine profiles and their control cohorts.

Permutational analysis of variance (PERMANOVA) was conducted in R 4.1.0 to compare the CP1_LI09 and CP2_LI09 cohorts for their bacterial and fungal microbiome compositions and compare them with their control cohorts (i.e., PC and NC cohorts). The same strategy was performed for the same microbiome composition comparisons of CP1_LI10 and CP2_LI10 cohorts and their control cohorts.

Similarity percentage (SIMPER) analysis was performed to determine the bacterial and fungal microbiome similarities within the LI09 and LI10 cohorts with different cytokine profiles. The same analysis was used to determine the bacterial and fungal microbiome dissimilarities between CP1_LI09 and CP2_LI09 cohorts. The same strategy was used for the comparisons of the bacterial and fungal microbiome dissimilarities between CP1_LI10 and CP2_LI10 cohorts.

Linear discriminant analysis (LDA) effect size (LEfSe) was carried out to compare the bacterial and fungi microbiome of CP1_LI09 and CP2_LI09 cohorts, respectively, to determine the bacteria and fungi associated with each of the two cohorts. The same analysis was performed to determine the bacteria and fungi associated with CP1_LI10 and CP2_LI10 cohorts.

Spearman test was used to determine the individual correlations between the cytokines in the cytokine profiles and the microbes associated with CP1_LI09, CP2_LI09, CP1_LI10, and CP2_LI10 cohorts.

Distance-based redundancy analysis (db-RDA) was performed to determine the effect of bacteria and fungi associated with CP1_LI09, CP2_LI09, CP1_LI10, and CP2_LI10 cohorts on the corresponding cytokine profiles.

Spearman test was carried out to investigate the correlations between the bacteria and fungi associated with CP1_LI09 cohort. The same strategy was performed for the bacteria and fungi correlations in CP2_LI09, CP1_LI10, and CP2_LI10 cohorts.

The correlations between the bacteria and fungi in the intestinal microbiome networks of LI09 and LI10 cohorts with different cytokine profiles were determined by co-occurrence network inference (CoNet) analysis. Five metrics, i.e., Bray–Curtis, Spearman, Pearson, mutual information, and Kullback–Leibler dissimilarities, were used to calculate the ensemble inference in CP1_LI09, CP2_LI09, CP1_LI10, and CP2_LI10 cohorts. The top 10 microbes (i.e., bacteria and fungi) with the largest numbers of correlations in the microbiome networks of the four cohorts were determined.

Network fragmentation calculations and generation of a null distribution were carried out in R as previously described (16) to explore the gatekeeper(s) in the microbiome networks of LI09 and LI10 cohorts with different cytokine profiles. Statistical significance was defined as the number of times a fragmentation score over that resulting from the removal of the bacteria or fungi observed within the null distribution.

Eleven serum cytokines were determined with significant differences among the LI09, LI10, PC, and NC cohorts, i.e., IL-1α, IL-2, IL-4, IL-5, IL-6, IL-12p70, IL-17A, M-CSF, MCP-1, MIP-3α, and RANTES (all p < 0.03), and they were used for the subsequent clustering analysis of cytokine profiles. Higher levels of IL-1α and M-CSF were determined in both LI09 and LI10 cohorts than in PC cohort, and IL-2 and IL-17A were greater in the LI10 cohort than in the PC cohort ( Supplementary Figure S1 ), suggesting IL-1α, M-CSF, IL-2, and IL-17A were enhanced by LI09 and/or LI10. In contrast, MCP-1, IL-5, and MIP-3α were determined with lower levels in both LI09 and LI10 cohorts than in the PC cohort ( Supplementary Figure S1 ), suggesting that they were suppressed by LI09 and LI10.

Two was determined as the most optimal number for clustering the cytokine profiles ( Figure 1A ), and two distinct cytokine profiles (i.e., CP1 and CP2) were determined in the LI09, LI10, PC, and NC cohorts ( Figure 1B ). All rats in the PC cohort were determined with CP1, and all rats in the NC cohort had CP2, suggesting that CP2 represented better immune status and was the “better cytokine profile.” Twenty-one rats in the LI09 cohort and 18 rats in the LI10 cohort were determined with CP1, while 19 rats in the LI09 cohort and 21 rats in the LI10 cohort had CP2 ( Figure 1B ).

Six out of the seven measured liver function variables, i.e., ALT, AST, ALP, TBA, TB, and GGT, were greater in the CP1_LI09 cohort than in the CP2_LI09 cohort ( Figure 2 ). By contrast, ALB was at similar level between the two cohorts.

Similarly, ALT, AST, ALP, TBA, TB, and GGT were determined at higher levels in the CP1_LI10 cohort than in the CP2_LI10 cohort ( Figure 3 ), while ALB was lower in the CP1_LI10 cohort than in the CP2_LI10 cohort ( Figure 3 ).

Ishak scoring system was used to help evaluate the liver histopathology, and in this study, a higher Ishak score represented greater liver injury severity. The Ishak score was greater in the CP1_LI09 cohort than in the CP2_LI09 cohort ( Figure 4A ), and the same score was greater in the CP1_LI10 cohort than in the CP2_LI10 cohort ( Figure 4B ). The rats with CP1 in LI09 and LI10 cohorts had more severe liver injury than those with CP2, further suggesting that CP2 was the “better cytokine profile” compared with CP1.

Firmicutes, Bacteroidota, and Verrucomicrobiota were determined as the three most abundant bacterial phyla in the LI09 and LI10 cohorts with different cytokine profiles.

At the family level, Lachnospiraceae and Bacteroidaceae were determined with the largest abundances in the bacterial microbiome of all the LI09 and LI10 cohorts with different cytokine profiles ( Figure 5 ). Akkermansiaceae was the third most abundant bacterial family in CP1_LI09, CP1_LI10, and CP2_LI09 cohorts, while Tannerellaceae was determined with the third most abundance in CP2_LI10 cohort ( Figure 5 ).

Shannon and Pielou indices were both similar in bacterial microbiome of the CP1_LI09, CP2_LI09, PC, and NC cohorts (both p > 0.09), while a significant difference was determined in observed species of the four cohorts (p < 0.001). CP1_LI09 and CP2_LI09 cohorts were found with similar observed species, while they were both greater than in the NC cohort ( Supplementary Figure S2A ). Likewise, Shannon and Pielou indices were both similar between CP1_LI10, CP2_LI10, PC, and NC cohorts (both p > 0.05), while observed species was different between the four cohorts (p < 0.001). No difference was found in the observed species of CP1_LI10 and CP2_LI10 cohorts, but they were both greater than in the PC and NC cohorts ( Supplementary Figure S2B ).

PERMANOVA analysis suggested that the bacterial composition was similar between CP1_LI09 and CP2_LI09 cohorts (R² = 0.028, p > 0.26), and they were both different from the PC and NC cohorts (p < 0.002) The same analysis determined a significant difference in the bacterial composition between CP1_LI10 and CP2_LI10 (R² = 0.063, p< 0.001), and they were both different from the PC and NC cohorts (p < 0.001).

SIMPER analysis determined that the similarity of bacterial phylotype abundances within CP1_LI09 was higher than that of CP2_LI09, i.e., 51% versus 46%. The same analysis determined a dissimilarity of 52% between CP1_LI09 and CP2_LI09 cohorts. Likewise, SIMPER analysis determined that the similarity of bacterial phylotype abundances was greater within CP1_LI10 (i.e., 51%) than within CP2_LI10 (i.e., 44%). The same analysis determined a dissimilarity of 55% between CP1_LI09 and CP2_LI09 cohorts.

LEfSe analysis determined five bacteria associated with CP1_LI09 cohort and one bacterium (i.e., Flavonifractor) associated with CP2_LI09 cohort ( Figure 6A ), among which Lachnospiraceae_UCG_006 was most associated with CP1_LI09 cohort. The same analysis revealed that 30 bacteria were associated with CP1_LI10 and CP2_LI10 cohorts ( Figure 6B ), among which Lachnospiraceae_NK4A136_group and Parabacteroides were most associated with CP1_LI10 and CP2_LI10 cohorts, respectively.

Ascomycota and Basidiomycota were determined as the two most abundant fungal phyla in all the LI09 and LI10 cohorts with different cytokine profiles.

At the family level, Aspergillaceae and Trichocomaceae were determined with most abundances in the LI09 and LI10 cohorts with different cytokine profiles ( Figure 7 ). Debaryomycetaceae, Hypocreales_Incertae_sedis, Trichosphaeriaceae, and Mrakiaceae were determined as the third largest abundances in the mycobiome of CP1_LI09, CP1_LI10 and CP2_LI09 and CP2_LI10 cohorts, respectively ( Figure 7 ).

Shannon and Pielou indices were both similar in the mycobiome between the CP1_LI09, CP2_LI09, PC, and NC cohorts (both p > 0.68), while observed fungal species was different in the four cohorts (p < 0.009). No difference was determined in the fungal observed species of CP1_LI09 and CP2_LI09, but they were both less than that in the NC cohort ( Supplementary Figure S3 ). Likewise, the three alpha diversity indices were all similar in the mycobiome of CP1_LI10, CP2_LI10, PC, and NC cohorts (all p > 0.53).

PERMANOVA analysis revealed that the mycobiome composition was similar between CP1_LI09 and CP2_LI09 cohorts (R² = 0.03, p > 0.22), but they were both different from PC and NC cohorts (p < 0.004). The same analysis determined a significant difference in the mycobiome composition between CP1_LI10 and CP2_LI10 cohorts (R² = 0.049, p < 0.007), and they were both different from PC and NC cohorts (p < 0.025).

SIMPER analysis determined that the similarity of fungal phylotype abundances within the CP1_LI09 cohort (i.e., 29%) was lower than that within the CP2_LI09 cohort (i.e., 34%). The same analysis determined a dissimilarity of 69% between CP1_LI09 and CP2_LI09 cohorts. Likewise, SIMPER analysis determined that the similarity of fungal phylotype abundances was lower within the CP1_LI10 cohort than that within the CP2_LI10 cohort, i.e., 34% versus 40%. The same analysis determined a dissimilarity of 65% between CP1_LI09 and CP2_LI09 cohorts.

LEfSe analysis determined seven fungi associated with CP1_LI09 and CP2_LI09 cohorts ( Figure 8A ), among which Meyerozyma and Penicillium were most associated with CP1_LI09 and CP2_LI09 cohorts, respectively. The same analysis determined that 12 fungi were associated with CP1_LI10 and CP2_LI10 cohorts ( Figure 8B ), among which Talaromyces and Aspergillus were associated with CP1_LI10 and CP2_LI10 cohorts, respectively.

Multiple correlations were determined between the cytokines in cytokine profile and the microbes associated with each of the CP1_LI09, CP1_LI10, and CP2_LI10 cohorts, except the CP2_LI09 cohort ( Figure 9 ). ASF356 and Meyerozyma were determined with more correlations in the CP1_LI09 cohort ( Figure 9A ). Lachnospiraceae_NK4A136_group, Meyerozyma, Stemphylium, and Talaromyces were determined with more correlations with the cytokines in the CP1_LI10 cohort, and Talaromyces seemed to have an opposite effect on IL-α and M_CSF when comparing with Meyerozyma and Stemphylium ( Figure 9C ). GCA_900066575 was negatively correlated with most cytokines in the cytokine profile of CP2_LI10 cohort ( Figure 9D ).

Some of the bacteria and fungi closely associated with CP1_LI09, CP2_LI09, CP1_LI10, and CP2_LI10 cohorts were determined to influence the cytokine profiles in the corresponding cohorts ( Figure 10 ).

Different correlations were determined between the bacteria and fungi associated with CP1_LI09, CP2_LI09, CP1_LI10, or CP2_LI10 cohorts. Three correlations were determined in the bacterial and fungi associated with the CP1_LI09 cohort, i.e., a positive correlation between Prevotellaceae_UCG_001 and Meyerozyma, negative correlations between Lachnospiraceae_UCG_006 and Alternaria and between ASF356 and Issatchenkia. By contrast, no correlation was determined between the bacteria and fungi associated with the CP2_LI09 cohort.

All the nine correlations were determined to be positive between six bacterial and six fungi associated with the CP1_LI10 cohort ( Figure 11 ). Defluviitaleaceae_UCG_011, Candidatus_Christensenellaceae, and Pygmaiobacter were positively correlated with Oidiodendron in the CP2_LI10 cohort.

CoNet results revealed 10 microbes with most correlations in the intestinal microbiome networks of LI09 and LI10 cohorts with different cytokine profiles ( Tables 1 , 2 ). Seven out of the top 10 microbes with most correlations in the CP2_LI09 cohort were not determined in the top 10 microbes in the CP1_LI09 cohort ( Table 1 ). Likewise, the top eight microbes with most correlations in the CP1_LI10 cohort were all not determined in the top 10 microbes in the CP2_LI10 cohort ( Table 2 ). Eubacterium and Clostridium were both determined with many correlations in the microbiome networks of both CP2_LI09 and CP2_LI10 cohorts but not in the probiotics-treated cohorts with CP1.

Among these bacteria and fungi, Hungatella, Papillibacter, and Myceliophthora were determined as gatekeepers in the network of the CP1_LI09 cohort (fragmentation analyses, all p < 0.05). In the CP2_LI09 cohort, Corynebacterium, Eubacterium, Papillibacter, and Cephalotrichiella were identified as the microbiome network gatekeepers (fragmentation analyses, all p < 0.05). By contrast, Udeniomyces and Pseudoflavonifractor were determined as the only gatekeepers in the microbiome networks of CP1_LI10 and CP2_LI10 cohorts, respectively (fragmentation analyses, both p < 0.03).

Fragmentation analysis showed that the fragmentation level of the microbiome network of CP1_LI09 cohort was lower than that of the CP2_LI09 cohort, i.e., 0.497 versus 0.575. Similarly, the fragmentation level was lower in the CP1_LI10 cohort (i.e., 0.416) than that in the CP2_LI10 cohort (i.e., 0.527).