A total of 64 healthy commercial crossbred boars (Duroc×Landrace×Yorkshire; 90 days 50.06 ± 3.75 kg, animals that had not received any antibiotics) were enrolled in a 60-day feeding trial. Animals were stratified by initial body weight and randomly allocated to two treatment groups, each comprising four pens with eight pigs per pen. The control group (CON) received the basal diet in dry pellet form, whereas the LFF group was provided fermented liquid feed prepared by mixing the basal diet with water (1:3, w/v) followed by daily batch fermentation. Diets were formulated according to NRC (2012) guidelines (Table 1). Feed was delivered using an automated centralized feeding system (Osborne, Shanghai, China) at 12:00 and 16:00, and pigs had unrestricted access to feed and water throughout the trial.

At the conclusion of the experiment, six pigs per treatment (twelve total) were randomly selected for sampling. Pigs were fasted for 8 h and deprived of water for 2 h prior to slaughter, then euthanized by electrical stunning followed by exsanguination. Blood (10 mL) was collected, allowed to clot at 4 °C, and centrifuged at 2500 r/min for 15 min to obtain serum, which was stored at −20 °C. Spleen tissues were excised immediately after slaughter, trimmed into ~1 cm³ sections on ice, placed into sterile cryovials, and snap-frozen in liquid nitrogen. To reduce intra-organ variability, samples used for sequencing were prepared as homogenized composites derived from both the left and right regions of the spleen (Three spleen samples were collected from each pig).

Dry feed was mixed with the bacterial inoculum at a solid-to-liquid ratio of 1:3 to prepare the fermentation substrate. The inoculum consisted of a compound microbial preparation dominated by Lactobacillus plantarum and supplemented with Saccharomyces cerevisiae. This microbial consortium possesses well-characterized probiotic properties, including the ability to inhibit common pathogenic bacteria in pigs, and is widely used to enhance microbial inoculation efficiency during the preparation of fermented liquid feed. Prior to inoculation, the microbial agent was dissolved or activated in a small volume of warm water (30–35 °C) and then evenly sprayed onto or incorporated into the feed–water mixture. Mechanical stirring was used to ensure uniform distribution of the inoculum. The final concentration of lactic acid bacteria at inoculation was 1.0 × 10⁷ CFU/mL.

The mixture was transferred into fermentation containers (Shanghai Bailun Biological Technology Co., Ltd., Shanghai, China), which were filled to no more than 70% of their total volume to allow adequate headspace. Excess air was expelled prior to sealing to create a semi-anaerobic environment. Fermentation was carried out at 26 ± 0.5 °C for 24–72 h. During fermentation, samples were collected every 12–24 h from multiple depths and locations using sterile pipettes or sampling spoons to measure pH, odor, microbial counts, and the presence of mold contamination.

The serum samples were thawed at 4 °C prior to analysis. The concentrations of interleukin-6 (IL-6), tumor necrosis factor-α (TNF-α), complement C3, immunoglobulin A (IgA), and interferon-γ (IFN-γ) were measured using commercial ELISA kits (Beyotime, China) according to the manufacturer’s instructions. Each sample was assayed in triplicate, and absorbance was recorded at 450 nm using a microplate reader. Cytokine concentrations were calculated based on the standard curves generated for each assay.

Total RNA was extracted using the TRIzol reagent (Invitrogen, CA, USA) according to the manufacturer’s protocol. RNA purity and quantification were evaluated using the NanoDrop 2000 spectrophotometer (Thermo Scientific, USA). RNA integrity was assessed using the Agilent 2,100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Then the libraries were constructed using VAHTS Universal V6 RNA-seq Library Prep Kit according to the manufacturer’s instructions. The transcriptome sequencing and analysis were conducted by OE Biotech Co., Ltd. (Shanghai, China).

The libraries were sequenced on an llumina Novaseq 6,000 platform and 150 bp paired-end reads were generated. About 47.9 M raw reads for each sample were generated. Raw reads of fastq format were firstly processed using fastp1 and the low quality reads were removed to obtain the clean reads. Then about 47.02 M clean reads for each sample were retained for subsequent analyses. The clean reads were mapped to the reference genome (https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/003/025/GCF_000003025.6_Sscrofa11.1/, Version Number: NCBI _ Sscrofa11.1) using HISAT22. FPKM3 of each gene was calculated and the read counts of each gene were obtained by HTSeq-count4. PCA analysis were performed using R (v 3.2.0) to evaluate the biological duplication of samples.

Differential expression analysis was performed using the DESeq2. p value < 0.05 and |log2foldchange| > 1 was set as the threshold for significantly differential expression gene (DEGs). Hierarchical cluster analysis of DEGs was performed using R (v 3.2.0) to demonstrate the expression pattern of genes in different groups and samples. The radar map of top 30 genes was drew to show the expression of up-regulated or down-regulated DEGs using R packet ggradar. Based on the hypergeometric distribution, GO (http://geneontology.org/), KEGG pathway (http://www.genome.jp/kegg/), Reactome and WikiPathways enrichment analysis of DEGs were performed to screen the significant enriched term using R (v 3.2.0), respectively. R (v 3.2.0) was used to draw the column diagram, the chord diagram and bubble diagram of the significant enrichment term.

A total of 30 mg of each sample was transferred to a 1.5mL EP tube. After adding 600 μL extract (containing L-2-chlorophenylalanine, 4 μg/mL), the sample was vortexed for 30 s, ultrasonically treated in an ice water bath for 10 min, and incubated for 1 h to precipitate the protein. The sample was then centrifuged at 12,000 r / min for 15 min. The supernatant was transferred to a new glass vial for analysis. Quality control (QC) samples were prepared by mixing equal amounts of supernatants from all samples.

LC–MS / MS analysis was performed using ACQUITY UPLC I-Class plus ultra-high performance liquid chromatography tandem QE high-resolution mass spectrometer: Chromatographic column: ACQUITY UPLC HSS T3 (100 mm × 2.1 mm, 1.8 um); column temperature: 45 °C; mobile phase: A-water (containing 0.1% formic acid), B-acetonitrile (containing 0.1% formic acid); flow rate: 0.35 mL / min; injection volume: 2 μL. Mass spectrometry conditions: ion source: ESI; sample mass spectrometry signal acquisition was performed in positive and negative ion scanning modes, respectively.

Raw metabolomics data were processed using Progenesis QI v3.0 (Nonlinear Dynamics, Newcastle, UK). Preprocessing steps included baseline correction, peak detection, peak integration, retention time alignment, and data normalization. Mass-matching tolerances were set to 5 ppm for precursor ions when referencing HMDB and Lipidmaps, and 10 ppm for the LuMet-Animal and METLIN databases; fragment ion tolerances were set to 10 ppm and 20 ppm, respectively. Metabolite annotation integrated multiple orthogonal features, including retention time (RT), accurate mass, MS/MS fragmentation patterns, and isotope distributions, and was performed against the Human Metabolome Database (HMDB), Lipidmaps (v2.3), METLIN, and the LuMet-Animal 3.0 local database.

Processed data were subsequently subjected to quality control procedures, including missing-value filtering, zero-value replacement, score-based filtering, and data merging. Ion features with mean missing (or zero) values exceeding 50% within any treatment group were removed. Remaining zero values were imputed with half of the minimum non-zero signal intensity across all samples (non-imputed values are documented in the accompanying missing-value matrix). Candidate metabolites were further filtered according to annotation scores generated by Progenesis QI. Annotation confidence was categorized into four levels (Level 1: RT ± 0.3 min and fragmentation score ≥ 45; Level 2: RT ± 0.3 min and fragmentation score < 45; Level 3: fragmentation score ≥ 45; Level 4: fragmentation score < 45). A threshold score of 36 (out of 80) was applied, and metabolites below this cutoff were excluded as unreliable identifications.

Multivariate statistical analyses were conducted using SIMCA 14.1 (Umetrics, Sweden). Principal component analysis (PCA) was first employed to examine sample clustering, distribution homogeneity, and potential outliers. Partial least squares discriminant analysis (PLS-DA) and orthogonal PLS-DA (OPLS-DA) were subsequently performed to identify discriminatory metabolic features between groups. Model performance was evaluated using R²X, R²Y, and Q² statistics, and model robustness was verified through 200-fold permutation testing. Variable importance in projection (VIP) scores were used to estimate the contribution of each metabolite to the discriminant model. Differential metabolites were defined as those with VIP > 1.0 and a false discovery rate P-adj < 0.05. Pathway enrichment analyses were conducted using the KEGG database (https://www.genome.jp/kegg/) via the MetaboAnalyst platform, and pathway significance and impact values were determined (p < 0.05).

To validate the accuracy of the sequencing results, a subset of differentially expressed genes was selected for quantitative real-time PCR (qRT-PCR) analysis (Each sample was analyzed in three technical replicate wells for PCR detection.). Total RNA was reverse-transcribed into cDNA using a reverse transcription kit (TIANGEN) according to the manufacturer’s instructions. qRT-PCR was performed on a Bio-Rad CFX96 Real-Time Detection System (Bio-Rad, Hercules, CA, USA) using the SYBR Premix Ex Taq kit (TaKaRa). Each 20 μL reaction contained 9 μL of SYBR Premix, 2 μL of cDNA, 1 μL each of forward and reverse primers, and 7 μL of ddH₂O. The amplification protocol consisted of an initial denaturation at 95 °C for 5 s, followed by 40 cycles of 95 °C for 5 s and 60 °C for 30 s, with a final step of 65 °C for 5 s and 95 °C for 5 s. Relative gene expression levels were calculated using the 2^-ΔΔCt method with GAPDH as the internal control. Primer sequences are listed in Supplementary Table S1.

Transcriptome and metabolome joint analysis was performed using R 3.2 and Python 3.9. First, differentially expressed genes (DEGs) and differentially expressed metabolites (DEMs) underwent separate normalization processes. DEGs underwent log₂ transformation based on FPKM values and low-expression genes were filtered out; DEMs underwent log₁₀ transformation based on peak area and Pareto scaling to reduce the impact of outliers. The mixOmics package (v6.20) was used to construct a cross-omics association model based on sparse partial least squares (sPLS). Outputs included a cross-omics correlation matrix and a weighted correlation network. Subsequently, a gene–metabolite correlation network was constructed. Spearman correlation coefficients were calculated for all DEGs and DEMs, retaining significant strong correlations (P-adj< 0.05 and |R| > 0.7). Select DEGs and DEMs with|R| > 0.8 for network visualization using Cytoscape (v3.9.0). Key nodes were identified based on degree centrality and betweenness centrality. Combining joint pathway enrichment analysis, MetaboAnalyst v5.0 was employed for metabolite enrichment in KEGG pathways, while clusterProfiler (v4.6) performed KEGG enrichment for DEGs. Significant pathways (p < 0.05) from both analyses were then filtered.

Qpcr and growth performance data data were statistically analyzed using GraphPad Prism software (v9.5). The difference between the two groups was evaluated by Student’s t test, and the results were expressed as mean ± standard deviation (mean ± SD). When p < 0.05, the difference was considered statistically significant.

During the whole experimental period, feeding fermented liquid feed (LFF) significantly improved the growth performance of pigs. Compared with the control group (CON), pigs in the LFF group had higher body weight at the end of the experiment and showed a significant increase in average daily gain and FCR (Table 2). In addition, the results of blood immunological indexes showed that the levels of serum IL-6, TNF-α, C3, IgA and IFN-γ in LFF group were significantly higher than those in CON group (Figure 1A), suggesting that fermented liquid feed could promote the growth and development of pigs and modulate the immune function of pigs.

In this study, a total of 12 cDNA libraries (CON-1 to CON-6 and LFF-1 to LFF-6) were constructed from spleen samples of the CON and LFF groups and subjected to high-throughput sequencing on the Illumina platform. On average, about 47.88 million and 47.93 million raw reads were obtained in the CON and LFF groups, respectively, and the quality assessment showed that the average GC content of the samples in the two groups exceeded 50% and the Q30 value was about 98%, which indicated that the sequencing data were highly reliable. After removing low-quality and redundant sequences, an average of 46.97 million and 47.06 million clean reads were obtained in the CON and LFF groups, respectively. Subsequently, HISAT2 (v2.2.1, http://github.com/infphilo/hisat2) was utilized to align the clean reads to the pig reference genome, and the alignment rate ranged from 91.89–93.77% (Supplementary Table S2). In conclusion, the quality of RNA sequencing data obtained in this study is reliable and can be used for subsequent bioinformatics analysis.

Principal component analysis (PCA) was performed based on the spleen transcriptomic profiles obtained from the CON and LFF groups. The first two principal components, PC1 and PC2, explained 22.38% and 17.08% of the total variance, respectively, and clearly distinguished the CON group from the LFF group (Figures 2A,B). This indicates that fermented liquid feed exerts a notable influence on the spleen transcriptome of pigs.Gene expression levels were then normalized using fragments per kilobase of transcript per million mapped reads (FPKM). The overall gene expression patterns between the two groups were found to be largely similar, as shown by the FPKM distribution plot.Differentially expressed genes (DEGs) were subsequently identified using the DESeq2 package (https://bioconductor.org/packages/release/bioc/html/DESeq2.html). A total of 184 DEGs were detected based on the criteria of |log₂ fold change| > 1 and p < 0.05, including 104 upregulated and 79 downregulated genes in the LFF group (Figure 2C). To validate the RNA-seq results, 12 DEGs were randomly selected for RT-qPCR analysis. The expression patterns of these genes were highly consistent with the RNA-seq data (Figure 2D), confirming the reliability and accuracy of the transcriptomic analysis.

To further explore the functions of differentially expressed genes (DEGs), gene ontology (GO) enrichment analysis was performed on 187 DEGs. The results showed that these DEGs were annotated to a total of 866 GO entries, including 232 biological process (BP), 88 cellular component (CC), and 113 molecular function (MF) entries. A total of 63 entries reached the significant level at a set threshold of p < 0.05, and the top 20 significantly enriched GO entries are shown in Figure 3A. Specifically, there were 31 significantly enriched entries in the BP category, involving cellular oxidant detoxification, leukocyte chemotaxis, neuropeptide signaling pathway and regulation of cell adhesion; Thirteen entries were significantly enriched in the CC category, including hemoglobin complex, haptoglobin-hemoglobin complex, and sodium channel complex; 19 entries were significantly enriched in the MF category, mainly involving chemokine Of the 31 notable entries in the BP category, many are closely related to immunity, such as immune response, chemokine-mediated signaling pathway, neurotrophic pathway, and the immune response. Further analysis showed that genes such as CXCL2, CXCL8, OSM and SLA-5 were closely associated with the immune response (Figure 3A).

KEGG pathway annotation was performed for the identified differentially expressed genes (DEGs). A total of 150 KEGG pathways were annotated (Figure 3B), with 27 pathways significantly enriched at p < 0.05. For instance, pathways including Dopaminergic synapse, Phagosome, Drug metabolism-cytochrome P450, and Parathyroid hormone synthesis, secretion and action were significantly enriched. Multiple genes were enriched in immune-related pathway entries, including the IL-17 signaling pathway, cAMP signaling pathway, and Aldosterone-regulated sodium reabsorption. Further analysis indicated that genes such as CXCL2, CXCL8, FOS, and FOSB were enriched in these immune-related pathways, suggesting their potential involvement in the regulation of host immune responses. KEGG pathway annotation was performed for the identified differentially expressed genes (DEGs).

Immune factor detection results revealed significant differences in serum immune factor levels between the CON and LFF groups, suggesting that feeding liquid fermented feed may influence metabolic patterns in pigs. To investigate this, we conducted metabolomics analysis on spleen samples from both groups using a high-resolution LC-MS platform. A total of 1,750 and 2,335 metabolites were identified in anion and ion modes, respectively. OPLS-DA analysis revealed that the two principal components, PC1 and PC2, explained 27.7% and 13.3% of the total variance, respectively. In the score plot, samples from the CON group and LFF group exhibited distinctly separate distribution regions, indicating that the two groups possessed clearly distinguishable metabolomic profiles (Figure 4A). Heatmap analysis further highlighted differences in metabolite abundance patterns between groups (Figure 4B).Based on this, differential metabolites were identified by setting screening thresholds: VIP > 1, |log2Fold Change| > 1, and p < 0.05. Results showed that in the anion mode, 34 differential metabolites were identified, including 23 significantly upregulated and 11 significantly downregulated; while 70 differentially expressed metabolites were identified in the cationic mode, comprising 38 significantly upregulated and 32 significantly downregulated metabolites (Figures 4C,D). Collectively, these findings demonstrate significant differences in the metabolic profiles between the CON and LFF groups, providing metabolomic support for the hypothesis that fermented liquid feed exerts regulatory effects on porcine physiological functions.

To explore the biological processes potentially involved in differential expression metabolites (DEMs) in vivo, KEGG enrichment analysis was performed. Results showed that a total of 96 KEGG pathways were annotated, with 71 significantly enriched at the p < 0.05 threshold (Figure 5A). Significantly enriched pathways primarily involved Biosynthesis of unsaturated fatty acids, Insulin secretion, Linoleic acid metabolism, Glycerophospholipid metabolism, Aldosterone synthesis and secretion, and Regulation of lipolysis in adipocytes.Notably, multiple immune-related signaling pathways were also significantly enriched, including NF-kappa B signaling pathway, Th1 and Th2 cell differentiation, Th17 cell differentiation, Chemokine signaling pathway, T cell receptor signaling pathway, and B cell receptor signaling pathway. Further pathway impact analysis revealed significant effects on Arachidonic acid metabolism and Glycerophospholipid metabolism (Figure 5B), suggesting these metabolic pathways may play a crucial role in LFF’s regulation of immunity and metabolism.

To further elucidate the relationship between transcriptomes and metabolomes, a joint analysis of differentially expressed genes (DEGs) and metabolites (DEMs) in the CON and LFF groups was conducted, establishing a gene-metabolite correlation network. The analysis indicates that fermented liquid feed may induce immune-related responses by regulating transcriptional levels and metabolic processes within the body. This network comprises 48 differentially expressed genes and 61 differential expression patterns, with 116 pairs of correlation relationships. Among these, SLC6A20, IL1RAPL2, SLA-5, OSM, and FAM209B exhibited significant positive correlations with Sarcophytrol J, N-Arachidonoyl Alanine, Glyceryl 5-Hydroxydecanoate, Gabapentin, Photinus Luciferin, Cyclohexanone, 4,8,12,15,19 -Docosapentaenoic Acid. Conversely, A2M exhibited negative correlations with 5-Dihydro-4,5-Dimethyl-2-(2-Methylpropyl)Thiazole, Methyl Dihydrojasmonate, 9-Tridecynoic Acid, and SM(d18:1/15:0). This suggests that FLF may influence the metabolic patterns of the host (Figure 6).