In this study, 209 Duroc gilts were generously donated by a commercial pig farm located in Sichuan, China, with an average age of 83.62 ± 5.37 days and an average weight of 33.68 ± 4.23 kg. The Duroc gilts were randomly allocated to 15 pens, with 12–14 pigs per pen, and each pen was equipped with a state-of-the-art Feed Intake Recording Equipment (FIRE) (RLX-096, Osborne, United States). Prior to the onset of the experiment, the pigs were allowed to acclimatize to the FIRE equipment over a duration of 1 week. The measurement process lasted for a total of 93 days, and each pig was assigned a unique radio frequency identification tag on the left ear for the FIRE system to detect. Upon entering the FIRE, the system automatically recorded pertinent parameters, including the number of visits to the feeder, feeding time, feed intake, and body weight. The swine resided within a climate-controlled environment, the temperature modulating between a mild 18 and a warm 25 degrees Celsius, accompanied by a relative humidity oscillating from 60% to a slightly damp 80%, with a 12 h light/dark cycle (lighting from 7 am to 7 pm). Throughout the trial, all the pigs were fed with the same corn-soybean commercial diet, conforming to the Chinese standard GB/T 5915-2020, and were not administered antibiotics or drugs. The primary constituents of the feed comprise maize, soya bean meal, and bran. It boasts a composition of 15% crude protein, 1.6% crude fat, 5% crude fiber, and 7.5% crude ash, in addition to noteworthy amounts of lysine (0.9%), calcium (0.8%), phosphorus (0.6%), and a dash of salt at 0.3%. Adequate supplies of clean water were provided ad libitum during the entire study period, and regular veterinary inspections were conducted to ensure that only healthy animals were selected for the trial.

We performed data quality control on the information gathered by the FIRE system, which were evaluated based on certain parameters: daily feed intake ranging from 0.5 to 4.5 kg, the number of feed intakes per day of 2 to 20 times, and daily feeding time between 5 min to 2 h. Any data that fell outside of these established ranges were identified as outliers and subsequently removed. After quality control, the dataset had reduced from the initial 105,202 entries to 92,305 entries. Each day’s weight was calculated based on the average measured body weight. Using the data that passed the quality control test, we then evaluated each pig’s total time spent at the feeder per day (TPD), number of visits to the feeder per day (NVD), average feed intake per minutes (AFI), measurement initiation body weight (W1), measurement termination body weight (W2), and average daily feed intake (ADFI). At the point when the pigs’ average weight reached 115 kg, their backfat thickness (BFT) in the tenth rib was determined using real-time B-mode ultrasound (MyLab^™X7, ESAOTE, Genova, Italy).

The average daily gain (ADG), ADFI, FCR, and RFI were calculated by the following model (22, 23):

Where TFI refers to total feed intake. AMW means average metabolic weight.

The HRC and LRC groups were comprised of the ten pigs with the highest and lowest RFI scores, respectively. Twenty select pigs were then transported to a commercial slaughterhouse for slaughter, where they remained on-site overnight, and fasted but with ad libitum access to water prior to being slaughtered the following morning. The selected live pigs were humanely slaughtered in accordance with the Live Pig Slaughter Guidelines (GB/T 17236-2019), which are approved by the General Administration of Quality Supervision, Inspection, and Quarantine of the People’s Republic of China and the Standardization Administration of the People’s Republic of China. After being stunned with carbon dioxide, they were bled for slaughter. The cecum was opened with sterile ophthalmic scissors, and the cecal contents were carefully collected with sterile forceps, immediately transferred to a sterile 50 mL centrifuge tube, and preserved in liquid nitrogen. Separate utensils were employed for each sample to ensure that there was no cross-contamination.

The QiAamp DNA Stool Mini Kit (Qiagen, Hilden, Germany) was employed to extract DNA from the cecal contents in line with the manufacturer’s instructions. Subsequently, approximately 30 ng of genomic DNA samples were taken and matching fusion primers were employed to establish an amplification PCR reaction system. The PCR amplification products were then purified with Agencourt AMPure XP magnetic beads (Beckman Coulter, Brea, CA), dissolved in Elution Buffer, labelled, and subjected to library construction. The fragment range and concentration of the library were determined using the Agilent 2,100 Bioanalyzer (Agilent, Santa Clara, CA, United States). Using the HiSeq platform, properly built libraries were sequenced based on the insert size. For sequence assembly, we used the Fast Length Adjustment of Short reads (FLASH) splicing program (v.1.2.11) to join the paired reads generated from paired-end sequencing into one sequence using overlapping relationships to produce tags of hypervariable regions (24). The splicing conditions were set such that the minimum match length was 15 bp, and the allowable mismatch rate in the overlapping region was 0. Afterwards, the spliced Tags were clustered into Operational Taxonomic Units (OTU s) (sequence similarity >97%) utilizing USEARCH (v.7.0.1090) software (25).

Calculation of the Alpha diversity analysis (including Chao1 index, Ace index, Shannon index, Simpson index and the Good’s coverage) was performed using MOTHUR (v.1.30.1) (26). Beta-diversity was quantified through the utilization of the Bray–Curtis distance metric, which was computed as resemblances (ANOSIM). Furthermore, the visualization of the results was accomplished by employing the QIIME2 package in the R software, employing principal coordinates analysis (PCoA) (27). Linear discriminant analysis effect size (LEfSe) analysis was performed using LEfSe software (28). PICRUSt2 (v.2.2.0b) was used for functional annotation of microbial communities (29). The rest of the graphs were drawn using the R (v 4.1.0). Statistically significant variations in the abundance of microbiota were observed at the phylum and genus levels, distinguishing between the HRC and LRC groups. This discrimination was determined by employing Student’s t-test and controlling for false discovery rate (FDR). Based on the 16S sequencing data, functional predictions based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database were performed using the PICRUSt2 software (30).

In this experiment, approximately 60 mg of cecal contents were extracted from each of the cryopreserved HRC and LRC groups. Ground the samples into homogenate using a tissue grinder. Then, added 400 mL of extraction solvent (methanol: water = 4:1) for extraction. Collected the supernatant by centrifugation (13,000 g, 15 min, 4°C) for further analysis. To separate and detect metabolites, a Waters 2D UPLC system (provided by Waters Corporation, United States) in combination with a tandem Q Exactive high-resolution mass spectrometer (manufactured by Thermo Fisher Scientific, MA, United States) was employed. The following chromatographic conditions were applied: employment of a BEH C18 column (100 mm length, 2.1 mm internal diameter, and 1.7 μm particle size; Waters, Milford, United States); the mobile phase for positive ion mode entailed an aqueous solution containing 0.1% formic acid (Solution A) and 100% methanol comprising 0.1% formic acid (Solution B). Conversely, the mobile phase for negative ion mode utilized an aqueous solution harboring 10 mM ammonium formate (Solution A) and 95% methanol containing 10 mM ammonium formate (Solution B). The elution was executed along the following gradient: 0–1 min, 2% Solution B; 1–9 min, 2–98% Solution B; 9–12 min, 98% Solution B; 12–12.1 min, a transition from 98% Solution B to 2% Solution B; 12.1–15 min, 2% Solution B. The sample injection volume was set at 5 μL, and the flow rate reached 0.36 mL/min. The column temperature was maintained at 45°C. Moreover, both positive and negative ion scan modes were employed to capture the signal emanating from the mass spectrometer analysis. Regarding the mass spectrometry settings, a stepped normalized collision energy protocol was utilized, employing values of 20, 40, and 60 eV respectively, while the spray voltage was set to 3.80 KV for positive mode and 3.20 KV for negative mode. Additional settings encompassed a capillary temperature of 320°C, an auxiliary gas heater temperature of 350°C, a sheath gas flow rate of 40 mL/min, and an auxiliary gas flow rate of 10 mL/min.

The raw data were preprocessed by metaX software for data preprocessing, statistical analysis, metabolite classification annotation, and functional annotation (31). The analysis was conducted utilizing the free online platform MetaboAnalyst 5.0¹ (32). Calculate the Variable importance in the projection (VIP) value from Partial Least Squares-Based Discriminant Analysis (PLS-DA). VIP ≥ 1.0, absolute fold change (FC) ≥ 2.0, p < 0.05 (Student’s t-test) were used as criteria for the selection of differential metabolites. The rest of the graphs were drawn using the R software (v.4.1.0). Pearson correlation analysis was used to determine the relationship between microbial communities and metabolites using the heatmap package (v.1.0.12) in R software.

The prediction of the heat map relied on a genome-scale metabolic model (GEM)-derived prognostic model. This prognostic model was built upon a logical regression model, trained on a meticulously crafted, GEM. It holed the capacity to anticipate the potential synthesis of each metabolite across various categorical tiers.

The data regarding growth performance and feeding behavior measurements were expressed as the mean ± standard deviation (SD). Statistical comparisons between groups were conducted using Student’s t-test with the aid of SPSS (v.22.0). For all results, the P criterion was as follows: *p < 0.05 was determined to be statistically significant, **p < 0.01 indicated a high level of statistical significance, and ***p < 0.001 denoted extremely significant.

First, we conducted statistical analysis on the production performance of the selected pigs (Table 1). The significant differences in RFI phenotypes prove the success of our grouping (p < 0.05). The results revealed that the LRC group exhibited lower FCR, ADFI (p < 0.001), and lower thicker backfat (p < 0.05) compared to HRC group.

Additionally, we assessed the feeding behavior of both groups, primarily focusing on the following metrics: TPD, NVD, and AFI. Interestingly, the LRC group exhibited markedly elevated levels of TPD (p < 0.05), contrasted by considerably and highly significantly reduced measures of AFI and NVD (p < 0.001) when juxtaposed with the HRC group.

A total of 1,111,597 readings were obtained from twenty specimens, averaging 55,579 readings per specimen. Following a thorough quality control and screening process, 1,288 OTUs were identified.

We have demonstrated the disparity in alpha diversity between the two groups, examining two facets: calculating the indices of bacterial richness estimators (Chao1, Ace) and computing the indices of bacterial diversity (Shannon, Simpson). The corresponding outcomes have been elaborated in Table 2. Notably, when comparing the LRC group to the HRC group, the flora’s abundance exhibited a significant difference (p < 0.01), whereas disparities in flora diversity were not statistically significant. Furthermore, the Good’s coverage analysis revealed a comprehensive coverage exceeding 99%, thus affirming the reliability of the findings.

In the beta diversity difference results, PCoA revealed distinct structural differences in the composition of the intestinal microbiota between the LRC and HRC group (p < 0.05) (Figure 1A). In light of the OTU analysis results derived from each sample, the widely utilized weighted unifrac index will be implemented to measure the dissimilarity coefficient between two specimens. As depicted in Figure 1B, a pronounced separation had been observed amongst the samples from both groups.

In the LRC group, three phyla, namely Bacteroidetes (22.53%), Firmicutes (51.28%), and Proteobacteria (19.19%), were dominant. Conversely, in the HRC group, Bacteroidetes (41.48%), Fusobacteria (32.45%), and Firmicutes (18.28%) played a pivotal role (Figure 2A). On a genus level analysis, the LRC group presented four prevailing genera: Escherichia (15.18%), Alloprevotella (7.69%), Eubacterium (6.78%), and Fusobacterium (5.26%). In comparison, the four dominant genera within the HRC group comprised of Fusobacterium (25.80%), Bacteroides (10.52%), Alloprevotella (5.51%), and Prevotella (4.64%) (Figure 2B).

As delineated in Table 3, we discerned six, nine, nine, and six markedly disparate microorganisms at the ranks of phylum, class, order, family, and genus correspondingly, each substantiating a FDR less than 0.05 (all differential microorganisms were shown in Supplementary Table S1).

Linear discriminant analysis scoring investigation has facilitated the detection of statistically substantial biomarkers amidst the groups, namely species exhibiting conspicuous dissimilarities across said groups. In the non-parametric Kruskal–Wallis rank sum test, the LEfSe analysis conducted at a threshold level of LDA ≥ 3.0 and p < 0.05 revealed 78 pathways that explained the distinguishing features between HRC and LRC pigs (Figure 3).

To explore into the functionalities of these distinctive microorganisms, we conducted KEGG pathway analyses using PICRUSt2. Results of considerable divergence were identified when mean relative abundances exceeded 1%, accompanied by P and FDR values below 0.05. A total of 18 KEGG pathways exhibited noteworthy dissimilarities between the LRC and HRC groups, as illustrated in Table 4.

Importantly, the LRC group epitomized elevated activity throughout pathways encompassing annamycin, pantothenate, CoA, and folate synthesis, as well as the formation of fatty acids, phenylalanine, tyrosine, and tryptophan. Furthermore, a marked amplification was noted in the metabolism of D-glutamine and D-glutamate, tryptophan, histidine, selenocompound, cysteine, and methionine.

Conversely, within the HRC group, there was a conspicuous augmentation in the relative abundance of streptomycin, lysine, valine, leucine, and isoleucine biosynthesis, protein transport mechanisms, the genesis of vancomycin group antibiotics, biotin metabolism, and the pentose phosphate pathway.

We performed a metabolomic analysis of cecal contents to understand how metabolites differed between the two groups. The 6,838 and 1,667 m/z features were detected in positive and negative ionization modes, respectively. We used the PLS-DA model to identify differential metabolites between the two groups. The parameters of the models: R² = 0.664 Q² = 0.382 for positive ionization mode; R² = 0.671 Q² = 0.426 for negative ionization mode (Figures 4A,D). The results showed that there was a clear separation between the two groups. Based on the criteria of p < 0.05, |Log₂FC| > 1, VIP > 1, we detected 390 (248 metabolites up and 142 down), 200 (97 metabolites up and 103 down) differential metabolites in positive and negative ionization modes (Figures 4B,E).

The top ten metabolites, which expressed prominent divergence in their levels, were elucidated in Table 5. Detailed metabolite information was listed in Supplementary Table S2. In addition, we did cluster analysis on the differential metabolite, and the findings showed good cluster quality (Figures 4C,F).

To delve deeper into the functionalities of these divergent metabolites, we conducted a KEGG enrichment analysis on the metabolites that exhibited noteworthy upregulation and downregulation in the LRC group in comparison to the HRC group. Surprisingly, the upregulated differentials were primarily concentrated within metabolic pathways, tryptophan metabolism, alpha-linolenic acid metabolism, linoleic acid metabolism, thiamine metabolism, and biosynthesis of unsaturated fatty acids (p < 0.05) (Figure 5). The down-regulated differential metabolites were mainly enriched in metabolic pathways, but not significantly (p > 0.05).

We compared the KEGG enrichment results of the two omics and found that they were both enriched in tryptophan metabolism. Specifically, six differential metabolites were identified as enriched in the tryptophan metabolic pathway. These metabolites include 4-(2-aminophenyl)-2,4-dioxobutanoic acid, 5-hydroxyindole-3-acetic acid, xanthurenic acid, 4-hydroxy-2-quinolinecarboxylic acid, serotonin, and 5-methoxyindoleacetic acid (as listed in Table 6). In order to explore deeper into the relationship between metabolites and microorganisms, we conducted a Spearman correlation analysis between the metabolites enriched in the tryptophan metabolic pathway (identified through the KEGG analysis) and the differential microbial genera (identified with LDA score > 4 through LEfSe analysis). The correlation analysis revealed noteworthy associations, as illustrated in Figure 6. Specifically, 4-(2-aminophenyl)-2,4-dioxobutanoic acid and 5-hydroxyindole-3-acetic acid both correlated positively with Megasphaera and Gemmiger respectively, as did 4-hydroxy-2-quinolinecarboxylic acid and xanthurenic acid with Megasphaera. Conversely, serotonin showed a pronounced negative correlation with Bacteroides, underscoring an inverse relationship between them (p < 0.001).