Integrated Omics Analysis of the Effects of Nano-Antimicrobial Peptides on the Intestinal Microbiota and Metabolome of Tibetan Sheep

1. Introduction

Tibetan sheep (Ovis aries), a unique and valuable livestock resource on the Qinghai-Tibet Plateau, are primarily distributed in the high-altitude pastoral regions of Qinghai, Tibet, and Gansu provinces, at elevations of around 4,000 meters above sea level. They play a vital role in plateau animal husbandry and the local economy [1]. Over long periods of natural selection, Tibetan sheep have developed a specialized intestinal microecological system that enables them to adapt to harsh environments, including extreme cold, hypoxia, and limited forage availability [2]. The energy acquisition, nutrient conversion, and immune homeostasis of Tibetan sheep heavily rely on a stable and functionally efficient intestinal microbiota.

The intestinal microbiota is a critical regulator of host health, influencing nutrient metabolism, immune function, and disease prevention [3,4,5]. Previous studies have shown that the microbiota of Tibetan sheep is primarily composed of dominant phyla such as Firmicutes, Bacteroidetes, Proteobacteria, and Actinobacteria [6,7]. These microorganisms play essential roles in energy metabolism, immune regulation, and the environmental adaptation of plateau ruminants [8,9,10,11]. Despite the importance of this microbial community, traditional breeding practices continue to rely heavily on antibiotics to regulate gastrointestinal health. Long-term use of antibiotics, however, can lead to bacterial resistance and food safety concerns, which challenges the sustainable development of animal husbandry [12].

Antimicrobial peptides (AMPs), natural host defense molecules, offer broad-spectrum antibacterial activity without the risk of inducing resistance, making them promising alternatives to antibiotics [13]. Paneth cells, a key source of AMPs, play a crucial role in regulating the intestinal microbiota and maintaining the mucosal barrier by secreting antimicrobial peptides such as α-defensins and lysozyme [5]. Their application has been proven effective in enhancing piglet weight gain, improving dairy cow milk quality, and optimizing ewe lactation performance [13,14,15]. Recent advancements in nanotechnology have further enhanced the stability and antibacterial efficacy of AMPs, making them even more effective in complex physiological environments and extending their in vivo half-life [16,17,18]. However, most existing studies focus on poultry and monogastric animals, with limited research on ruminants, particularly Tibetan sheep, and the mechanisms of intestinal microecological responses to nano-antimicrobial peptides (NAPs) remain unclear.

To address this gap, this study aimed to evaluate the effects of dietary NAPs on the intestinal microbiota of Tibetan sheep, using 16S rRNA gene sequencing and non-targeted metabolomics. The goal was to provide a theoretical foundation for green breeding of Tibetan sheep in alpine pastoral areas and the development of novel antibiotic-free feed additives. In this study, the rumen, small intestine, and rectum were selected as sampling sites. The small intestine is the primary site for nutrient digestion and absorption, and its mucosal structure and microbiota directly influence digestive efficiency and host health. As a specialized digestive organ in ruminants, the rumen plays a central role in feed degradation, volatile fatty acid production, and microbial metabolism. The rectum, located at the distal end of the gastrointestinal tract, serve as an indirect indicator of the microecological status of the distal gut. Together, these three sites represent distinct functional regions of the gastrointestinal tract in Tibetan sheep, allowing for a comprehensive understanding of the effects and mechanisms of nano-antimicrobial peptides in the gastrointestinal tract.

2. Materials and Methods

2.1. Ethics Approval

All experimental procedures were approved by the Institutional Animal Care and Use Committee (IACUC) of the College of Veterinary Medicine, Qinghai University (Approval No: SL-2023042; Approval Date: 2023.06.29). Every effort was made to minimize animal suffering and to reduce the number of animals used.

2.2. Animals and Experimental Design

Eighteen non-pregnant Tibetan sheep, weighing 35 ± 3 kg, were obtained from the experimental station of Qinghai University, Qinghai, China. The sheep were randomly assigned to two groups (n = 9 per group): Group A (Control group): fed a basal diet consisting of 500 g concentrate and 300 g hay per sheep/day; Group B (Experimental group): fed the same basal diet as Group A supplemented with 0.125 g of NAPs per sheep/day.

Sheep were housed in standardized individual pens with controlled ventilation, natural lighting, and a regulated environment. All animals had ad libitum access to clean drinking water. The facilities were cleaned and disinfected daily to maintain high hygiene standards and minimize environmental stress. Health status was monitored daily by a veterinarian throughout the 30-day experimental period.

2.3. Sample Collection

After 30 days of feeding, all sheep were slaughtered for sample collection. The animals were anesthetized via intravenous injection of sodium pentobarbital (20 mg/kg) and then exsanguinated without prior fasting in order to preserve the natural gastrointestinal contents.

Following laparotomy, the small intestine, rumen, and rectum were isolated. Contents from these three sites were collected and systematically labeled. For Group A, samples were designated as A-a (small intestine), A-b (rumen), and A-c (rectum); Group B samples were similarly labeled as B-a, B-b, and B-c. Each sample (approx. 2 mL) was immediately transferred into sterile 2.5 mL cryovials, snap-frozen in liquid nitrogen, and subsequently transported on dry ice to Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China) for 16S rRNA gene sequencing and metabolomic analysis.

2.4. Bacterial 16S rRNA Gene Sequencing and Analysis

Total genomic DNA was extracted from the gastrointestinal contents using the HiPure Stool DNA Kit (Meiji Biotechnology Co., Ltd., Guangzhou, China) according to the manufacturer’s instructions. The concentration and purity of the extracted DNA were assessed using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, USA) to ensure sufficient quality for downstream applications. The V3–V4 hypervariable region of the bacterial 16S rRNA gene was amplified using the universal primers 341F (5′-CCTACGGGNGGCWGCAG-3′) and 806R (5′-GACTACHVGGGTATCTAATCC-3′). PCR amplification was performed using a high-fidelity enzyme system in triplicate to ensure technical consistency. The resulting amplicons were purified using AMPure XP Beads (Beckman Coulter, USA) and quantified with a Qubit 3.0 Fluorometer (Thermo Fisher Scientific, USA). Sequencing libraries were constructed using the Nextera XT Index Kit (Illumina, USA) with unique sample barcodes. After library quality validation, 2 × 300 bp paired-end sequencing was performed on the Illumina MiSeq platform. Raw sequencing data were processed using QIIME 2 (version PE250). Raw paired-end reads were first filtered to remove low-quality sequences and then assembled into tags. These tags underwent further quality control to generate clean tags. Subsequently, sequences were clustered into OTUs, during which chimeric sequences were identified and removed to yield effective tags. Taxonomic classification was assigned to the representative sequence of each OTU using the SILVA database (Release 138) with a confidence threshold of 0.8. Based on the OTU abundance profile, subsequent analyses, including alpha and beta diversity, as well as functional prediction, were performed.

2.5. LC-MS Non-Targeted Metabolomics Detection and Analysis

Samples were thawed and extracted using a pre-cooled mixture of methanol, acetonitrile, and water (2:2:1, v/v). Following vortexing and ultrasonic extraction (30 min, 4 °C), the mixture was kept at −20 °C for 10 min. After centrifugation at 14,000 × g (20 min, 4 °C), the supernatant was collected and evaporated to dryness under vacuum. For final analysis, the residue was redissolved in 100 μL of 50% acetonitrile, vortexed, and centrifuged at 14,000 × g for 15 min at 4 °C to remove any remaining particulates. The final supernatant was collected for injection. Quality Control (QC) samples were prepared by pooling equal aliquots from each test sample to equilibrate the LC-MS system, monitor instrument performance, and ensure experimental reproducibility.

Metabolite separation was performed using an Agilent 1290 Infinity UHPLC system equipped with a HILIC column. The column temperature was maintained at 25 °C, with a flow rate of 0.5 mL/min and an injection volume of 2 μL. Mobile phase A consisted of water containing 25 mM ammonium acetate and 25 mM ammonium hydroxide, while mobile phase B was 100% acetonitrile. The gradient elution program was optimized as follows: 0–0.5 min, 95% B; 0.5–7 min, 95%–65% B; 7–8 min, 65%–40% B; 8–9 min, 40% B; 9–9.1 min, 40%–95% B; and 9.1–12 min, 95% B. Samples were kept in a 4 °C autosampler throughout the analysis, with QC samples inserted into the sequence at regular intervals.

Mass spectrometry was conducted using an AB Triple TOF 6600 system to acquire both primary (MS1) and secondary (MS2) spectra. Electrospray ionization (ESI) source parameters were set as follows: Gas1/Gas2, 60 psi; curtain gas (CUR), 30 psi; ion source temperature, 600 °C; and spray voltage, ±5500 V. The MS1 m/z range was 60–1000 Da with an accumulation time of 0.20 s/spectra. The MS2 m/z range was 25–1000 Da with an accumulation time of 0.05 s/spectra. Data-dependent acquisition (DDA) was performed using the Information-Dependent Acquisition (IDA) mode with a declustering potential of ±60 V and a collision energy of 35 ± 15 eV. Isotope ions were dynamically excluded within 4 Da, and 10 fragment spectra were collected per scan.

To facilitate subsequent analysis, raw spectral data were first converted to the .mzML standard using ProteoWizard. Then the XCMS algorithm was utilized to perform automated peak alignment and retention time adjustment, effectively normalizing the extracted peak areas across all samples..

2.6. Statistical Analysis

For alpha diversity indices, Welch’s t-test was used to evaluate differences between groups. Principal Coordinate Analysis (PCoA) based on Bray-Curtis distances was performed to visualize microbial community shifts, with significance assessed via PERMANOVA (P < 0.05). To identify specific microbial biomarkers, pairwise intergroup comparisons were conducted using the Wilcoxon rank-sum test. Linear Discriminant Analysis (LDA) Effect Size (LEfSe) was applied to estimate the effect size of each differentially abundant taxon, with an LDA score threshold of > 2.0 and P < 0.05. Functional metagenomes were predicted based on the OTU abundance table using PICRUSt2 software to infer Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment.

Metabolite annotation was performed by matching accurate mass and MS/MS fragments against public databases, including MassBank, METLIN, and MoNA, supplemented by an in-house secondary mass spectrometry database. Differentially abundant metabolites (DAMs) were identified based on the following criteria: Variable Importance in Projection (VIP) ≥ 1 from the OPLS-DA model and P < 0.05 from Student’s t-test, and |log₂FC| > 1 (fold change > 2 or < 0.5).

Metabolic pathway enrichment analysis was conducted via the KEGG database, with significance determined by a P value or false discovery rate (FDR) < 0.05 following -log10 transformation. To integrate the two datasets, an O2PLS-DA model was utilized to identify common variation dimensions between the microbiome and metabolome. Spearman correlation analysis was subsequently performed using the pheatmap package in R (version https://CRAN.R-project.org/package=pheatmap) to elucidate associations between differential microbial genera and metabolites.

3. Results

3.1. Analysis of 16S rRNA Gene Sequencing Results

To investigate the effects of NAPs supplementation on the gastrointestinal microbiota of Tibetan sheep, 16S rRNA gene sequencing was performed on gastrointestinal samples collected from Groups A and B after 30 days of feeding. After quality control and filtering, more than 80% of raw reads were retained as high-quality sequences for each sample. The number of effective tags exceeded 40,000 per sample, indicating sufficient sequencing depth for reliable downstream analyses. Statistical evaluation of Operational Taxonomic Unit counts (OTUs) and sequencing depth metrics was conducted (Figure 1). Variations in OTUs and tag counts were observed among samples, reflecting differences in microbial richness across groups. These results established the basis for subsequent analyses of microbial community structure and diversity.

3.2. Alpha Diversity Analysis

The OTUs in each sample were ranked according to their relative abundance to construct rank–abundance curves. Along the horizontal axis, the length of the curve reflects species richness, with a longer curve indicating a greater number of OTUs. Along the vertical axis, the slope of the curve represents species evenness; a flatter curve indicates a more even distribution of species abundance within the sample. As shown in Figure 2, samples A-a and B-a demonstrated relatively steeper slopes compared with other groups, indicating reduced species evenness. In contrast, the remaining groups exhibited flatter curves, suggesting a more homogeneous distribution of microbial taxa.

Based on Welch’s t-test, alpha diversity analysis using the Sob, Chao1, and ACE indices (Table 1) revealed significant differences (P < 0.05) in community richness between groups A-a and B-a. However, no significant differences in species richness were observed between groups A-b and B-b or between groups A-c and B-c.

3.3. Beta Diversity Analysis

PCoA was performed to evaluate differences in microbial community structure among groups. In the PCoA plot, samples positioned closer together indicate more similar microbial community compositions, whereas greater separation reflects increased compositional dissimilarity. As shown in Figure 3 and Table 2, a significant difference in microbial community structure was observed between groups A-a and B-a (P < 0.05). In contrast, no significant differences were detected between groups A-b and B-b or between groups A-c and B-c (P > 0.05). These results were consistent with the alpha diversity analysis.

3.4. Differential Microbial Composition Analysis

Based on the results of the alpha and beta diversity analyses, significant differences in microbial richness and community structure were observed between groups A-a and B-a, whereas no significant differences were detected between groups A-b and B-b or between groups A-c and B-c. Therefore, subsequent microbial composition analysis focused on groups A-a and B-a. The results are presented in Table 3. The top 10 dominant genera in both groups showed a largely similar taxonomic composition, with Candidatus_Saccharimonas, Christensenellaceae_R7_group, Ruminococcus, NK4A214_group, and Saccharofermentans constituting the core microbiota. However, notable differences were observed in their relative abundances. In group A-a, Candidatus_Saccharimonas was the most abundant genus, whereas in group B-a, Saccharofermentans became the predominant genus, and Candidatus_Saccharimonas decreased to second place. These findings suggest that NAPs supplementation did not alter the core taxonomic composition of the intestinal microbiota in Tibetan sheep but significantly reshaped the relative abundance structure of dominant genera.

Based on the significant differences observed in alpha and beta diversity between groups A-a and B-a, linear discriminant analysis effect size (LEfSe) was performed to identify taxa that were differentially enriched between the two groups. The results are shown in Figure 4. The LEfSe analysis revealed that Bacillota, Bacteria, Bacteroidota, and Pseudomonadota were significantly enriched in group B-a. In contrast, Patescibacteria, Archaea, Methanobacteriota, Myxococcota, and Nitrospirota were significantly enriched in group A-a.

3.5. Functional Prediction of Microbial Communities

Based on the significant differences observed in alpha diversity, beta diversity, and LEfSe analyses, functional prediction of the microbial communities was conducted. The results are shown in Figure 5. At KEGG, metabolism was the predominant functional category in both groups A-a (Figure 5A) and B-a (Figure 5B). Within this category, carbohydrate metabolism showed the highest relative abundance, followed by metabolism of cofactors and vitamins and amino acid metabolism. Genetic information processing represented the second most abundant functional module in both groups. Although the overall functional profiles were similar between A-a and B-a, differences were observed in the relative abundances of several metabolic pathways, suggesting that NAPs supplementation may have modulated functional potential without fundamentally altering core metabolic functions.

The above results indicate that although the overall taxonomic composition of the core microbiota remained stable between groups A-a and B-a, the relative abundances of dominant taxa were altered. This suggests that NAPspeptide supplementation may influence gastrointestinal metabolic potential by modulating the proportional distribution of key microbial taxa. Functional prediction analysis further demonstrated that no significant differences were observed in the overall functional gene profiles between the two groups. However, variations in the relative abundances of specific metabolic pathways were detected, implying potential shifts in microbial functional capacity. To further elucidate the impact of NAPs on gastrointestinal metabolism in Tibetan sheep, non-targeted metabolomic analysis was subsequently performed. The results are presented below.

3.6. Analysis of Non-Targeted Metabolomics Results

To evaluate data quality and analytical stability, the combined positive and negative ion datasets were subjected to principal component analysis (PCA), including QC samples (Figure 6). The PCA score plot showed that QC samples were tightly clustered, with overlapping distribution in the multivariate space, indicating good analytical reproducibility and system stability throughout the experiment. Metabolite identification and annotation were subsequently performed (Table 4). A total of 18,102 metabolite features were detected, including 9,793 in positive ion mode and 8,309 in negative ion mode. Among these, 3,916 positive ion features and 3,521 negative ion features were successfully annotated to known metabolites.

3.7. Differential Metabolite Profiling Between Groups

Based on the abundance matrix containing 7,437 annotated metabolite features retained across all samples, differential metabolite analysis was performed between Group A (control) and Group B (experimental). Metabolites were screened using the criteria of log₂ fold change (FC＞1) or (FC＜-1)and P < 0.05. The distribution of differential metabolites is presented as volcano plots in Figure 7. In the plots, red dots represent upregulated metabolites, whereas yellow dots indicate downregulated metabolites. Overall, more metabolites were downregulated than upregulated in all comparisons.. Notably, several downregulated metabolites exhibited high variable importance in projection (VIP) values, suggesting that these metabolites contributed substantially to group discrimination and may represent key metabolites distinguishing Group A from Group B.

A statistical summary of differentially expressed metabolites between groups is presented in Figure 8. In the comparison between A-a and B-a, 444 metabolites were upregulated and 1,424 were downregulated in group B-a relative to group A-a. Similarly, in the A-b vs. B-b comparison, 199 metabolites were upregulated and 1,114 were downregulated. In contrast, the A-c vs. B-c comparison showed 432 upregulated and 203 downregulated metabolites. Overall, the number of downregulated metabolites exceeded that of upregulated metabolites in the A-a vs. B-a and A-b vs. B-b comparisons, whereas the opposite trend was observed in the A-c vs. B-c comparison.

Differential metabolite analysis further identified key discriminative metabolites between groups A and B based on variable importance in projection (VIP) scores derived from the OPLS-DA model (Figure 9). In the A-a vs. B-a comparison, glycocholic acid, linoleic acid, glycerophosphocholine, oleic acid, and lysophosphatidylcholine (18:2) were enriched in group B-a, whereas cholic acid, 1-stearoyl-2-hydroxy-sn-glycero-3-phosphocholine, 1-palmitoyl-sn-glycero-3-phosphocholine, deoxycholic acid, D-desthiobiotin, trihydroxycholan-24-oic acid, and 15-cyclohexylpentanorprostaglandin F2α were enriched in group A-a. In the A-b vs. B-b comparison, hexadecanedioic acid, leptomycin B, and (−)-O-acetyl-D-mandelic acid were enriched in group B-b, whereas 15-cyclohexylpentanorprostaglandin F2α, pentadecanoic acid, enterolactone, heptadecanoic acid, 3α,7α-dihydroxy-12-oxocholanoic acid, and cholesterol sulfate were enriched in group A-b. Similarly, in the A-c vs. B-c comparison, (2E,4E)-12-hydroxy-13-(hydroxymethyl)-14-methoxy-3,5,7-trimethyl-14-oxotetradeca-2,4-dienoic acid, 15-ketoprostaglandin F2α, palmitoleoyl ethanolamide, leptomycin B, and aurantio-obtusin were enriched in group B-c, whereas 1H-pyrazole-3-carboxylic acid, pheophorbide a, antimycin A, and 1-palmitoyl-2-oleoyl-phosphatidylglycerol were enriched in group A-c.

3.8. KEGG Pathway Enrichment Analysis of Differential Metabolites Between Groups

Differential metabolites were mapped to the KEGG database for pathway enrichment analysis. A hypergeometric test was applied to identify pathways significantly enriched in differential metabolites relative to the overall metabolite background. As shown in Figure 10A, several pathways were significantly enriched in the A-a vs. B-a comparison, including choline metabolism in cancer, histidine metabolism, and primary bile acid biosynthesis. In contrast, no significantly enriched pathways were detected in the A-b vs. B-b comparison (Figure 10B). Only a limited number of pathways were enriched in the A-c vs. B-c comparison (Figure 10C), including alpha-linolenic acid metabolism and biosynthesis of unsaturated fatty acids.

Non-targeted metabolomics analysis indicated that, unlike the significant pathway enrichment observed in the A-a vs. B-a comparison, no significantly enriched pathways were detected between A-b and B-b. This suggests that the differential metabolites identified in these two subgroups did not converge into coordinated pathway-level alterations but were more likely to represent isolated metabolic fluctuations with limited systemic impact. In contrast, a small number of pathways were enriched in the A-c vs. B-c comparison, including alpha-linolenic acid metabolism and biosynthesis of unsaturated fatty acids. Combined with the profile of differential metabolites in this subgroup, these findings suggest potential localized modulation of lipid metabolism. However, no extensive pathway activation or inhibition was observed. Overall, these results indicate that the most pronounced and coordinated metabolic alterations occurred in the A-a vs. B-a comparison, highlighting this subgroup as the primary focus of metabolic regulation in the present study.

3.9. Integrated Analysis of 16S rRNA Gene Sequencing and Non-Targeted Metabolomics

To further investigate the overall associations between the intestinal microbiota (genus level) and the metabolomic profiles of Tibetan sheep following NAPs supplementation, an integrative multi-omics analysis was conducted. An orthogonal two-way partial least squares (O2PLS) model was constructed to evaluate the covariance structure between the microbiome and metabolome datasets. The contribution statistics of the O2PLS model are summarized in Table 5.

The shared variation between the microbiome and metabolome datasets was extracted using the O2PLS model. As shown in Figure 11A, genera such as Syntrophococcus and Family_XIII_UCG-001 were distributed on the left side of the loading plot, whereas Intestinimonas and Oscillospira were located on the right side, indicating distinct contributions of different genera to the integrative model. Similarly, in Figure 11B, metabolites such as prenosterone and ethylmorphine were positioned in the upper right quadrant, while talkoxydir and sulfate were located in the lower right quadrant, suggesting functional differentiation among metabolites within the integrated association structure. To further clarify the association patterns between the intestinal microbiota (genus level) and metabolomic profiles, a Spearman correlation heatmap was constructed using the R package pheatmap based on genus-level microbial abundance and differential metabolites. The results are shown in Figure 12. Metabolites and genera clustered within the same color blocks exhibited similar correlation patterns, implying potential involvement in related metabolic pathways or physiological processes. Terbutaline and 2-dhahba (dmed-fahfa) showed significant positive correlations (P < 0.05) with certain genera, including Lacrimispora. In contrast, most metabolites, such as prostaglandin D3 and chenodeoxycholic acid, displayed significant negative correlations (P < 0.05) with genera including Lacrimispora and Clostridium_methylpentosum.

The correlation network exhibited a clear modular structure, in which distinct microbial genera and metabolites formed relatively independent interaction clusters, reflecting the complexity and specificity of microbiome–metabolome associations in the intestinal ecosystem of Tibetan sheep. Within the network, genera such as Lachnospiraceae_NK4A136_group, Dialister, and Aerotrigona functioned as highly connected nodes, suggesting their potential roles as key hubs in host metabolic regulation. Lachnospiraceae_NK4A136_group is known to participate in short-chain fatty acid (SCFA) production, and its positive correlations with specific metabolites may indicate its involvement in dietary fiber fermentation and subsequent metabolic processes. Dialister, an acid-producing genus, showed associations with metabolites potentially related to intestinal acidification and barrier function. Similarly, the positive correlations observed between Aerotrigona and several metabolites may reflect its potential contribution to carbohydrate metabolism or host energy regulation.

Figure 13. Integrated microbiome–metabolome correlation network at the genus level. Note: Blue nodes represent metabolites, and orange nodes represent microbial genera. Red edges indicate significant positive correlations, whereas blue edges indicate significant negative correlations. Edge thickness reflects the strength of the correlation coefficient, and node size corresponds to the degree of connectivity (number of associations).

Figure 13. Integrated microbiome–metabolome correlation network at the genus level. Note: Blue nodes represent metabolites, and orange nodes represent microbial genera. Red edges indicate significant positive correlations, whereas blue edges indicate significant negative correlations. Edge thickness reflects the strength of the correlation coefficient, and node size corresponds to the degree of connectivity (number of associations).

4. Discussion

Based on 16S rRNA gene sequencing and non-targeted metabolomics, this study systematically evaluated the effects of NAPs supplementation on the microbial community structure and metabolic profiles of different gastrointestinal sites in Tibetan sheep and further explored the potential associations between the microbiome and metabolome. Quality control results demonstrated adequate sequencing depth and a high proportion of high-quality reads that met the technical specifications for gastrointestinal microbial research in ruminants [19], indicating robust and reliable data for downstream analyses. Rank–abundance and alpha diversity analyses revealed significant differences in microbial richness and evenness in the small intestinal contents between the experimental and control groups. Consistently, beta diversity analysis showed clear separation of small intestinal microbial communities between groups, whereas no significant differences were observed at other gastrointestinal sites. These findings suggest that NAPs supplementation exerts site-specific modulation of the intestinal microbiota in Tibetan sheep, with the small intestine appearing particularly responsive. Although the dominant core genera were largely conserved between groups, shifts in their relative abundances were observed. Functional prediction analysis further indicated that while the overall functional gene profiles were not significantly altered, variations in the relative abundance of specific metabolic pathways were detected. Together, these results imply that NAPs may influence gastrointestinal microbial function primarily through modulation of taxonomic structure rather than through large-scale restructuring of the core microbial community.

Metabolomic analysis further demonstrated distinct metabolic profiles between the experimental and control groups. The tight clustering of QC samples confirmed the stability and reproducibility of the analytical platform [20]. A total of 18,102 metabolite features were detected, and differential analysis was performed based on 7,437 annotated compounds consistently retained across all samples. Using |log₂FC| ≥ 1.5 and P < 0.05 as screening criteria, differential metabolites were primarily enriched in pathways related to bile acid metabolism, fatty acid metabolism, and phospholipid metabolism. In Group A, higher levels of cholic acid, deoxycholic acid, and trihydroxycholan-24-oic acid were observed. These bile acids are central components of the bile acid metabolic pathway and play essential roles in cholesterol conversion, lipid digestion, and nutrient absorption. Previous studies have also highlighted the close relationship between bile acid homeostasis, intestinal microecology, and lipid metabolism [21]. Additionally, elevated levels of phosphatidylcholine-related metabolites, including 1-stearoyl-2-hydroxy-sn-glycero-3-phosphocholine and 1-palmitoyl-sn-glycero-3-phosphocholine, may reflect maintenance of membrane structural integrity and lipid transport capacity [22]. The higher abundance of D-desthiobiotin and related biotin derivatives further suggests active participation in fundamental metabolic processes such as fatty acid synthesis and amino acid metabolism. Collectively, these findings indicate that the metabolic profile of Group A was characterized by relatively stable bile acid and phospholipid metabolism. In contrast, Group B exhibited increased levels of unsaturated fatty acids, including linoleic acid and oleic acid. These metabolites are closely associated with lipid metabolic balance, membrane fluidity, and inflammatory regulation [23]. Such changes may reflect adaptive adjustments in lipid metabolism following NAPs supplementation. Moreover, the elevated abundance of glycocholic acid and phospholipid-related intermediates, such as lysophosphatidylcholine (18:2) and glycerophosphocholine, suggests modulation of phospholipid turnover and remodeling processes, potentially influencing membrane dynamics and lipid-mediated signaling [24]. Taken together, the coordinated alterations in bile acid and phospholipid metabolism in Group B imply a restructuring of lipid metabolic networks rather than a complete disruption of metabolic homeostasis.

The differential metabolites were analyzed for pathway enrichment using a hypergeometric test. The results revealed significant enrichment in several pathways in the A-a vs. B-a comparison, such as choline metabolism in cancer, histidine metabolism, and primary bile acid biosynthesis. In contrast, no notable pathway enrichment was found in the A-b vs. B-b comparison, and only a few pathways, including alpha-linolenic acid metabolism and unsaturated fatty acid biosynthesis, were enriched in the A-c vs. B-c comparison. Interestingly, a clear functional relationship between the metabolites and enriched pathways was observed in the A-a vs. B-a comparison.. Group A-a exhibited higher levels of primary bile acids and phosphatidylcholine-related metabolites, which may reflect relatively stable bile acid biosynthesis and membrane-associated lipid metabolism. In contrast, Group B-a showed elevated levels of conjugated bile acids, phospholipid metabolism intermediates, and unsaturated fatty acids, suggesting a potential shift in metabolic flux within these pathways. These findings suggest that bile acid and lipid metabolic pathways may be the central nodes in the metabolic response to NAPs supplementation. Future studies integrating targeted metabolite quantification and functional validation experiments will be necessary to determine whether these metabolic alterations directly mediate the physiological effects observed.

From the perspective of microbiome–metabolome integration, the dynamic coupling between intestinal microbial communities and host metabolic networks is widely considered a key mechanism underlying physiological homeostasis [25]. To characterize the co-variation structure between these two omics layers, an O2PLS model combined with Spearman correlation analysis and heatmap visualization was applied to explore potential associations between microbial genera and metabolites. The results indicated that terbutaline and 2-dhahba (dmed-fahfa) were significantly positively correlated with the genus Lacrimispora. This association suggests that Lacrimispora may be linked to fatty acid– or steroid-related metabolic processes. As a member of the phylum Firmicutes, Lacrimispora is commonly regarded as an intestinal fermentative bacterium and has been reported to participate in fatty acid metabolism and related pathways [26]. Therefore, fluctuations in its abundance may be associated with changes in lipid-related metabolic activity. In contrast, several metabolites, including prostaglandin D3 and chenodeoxycholic acid, showed significant negative correlations with genera such as Lacrimispora and Clostridium_methylpentosum. Chenodeoxycholic acid, an important bile acid component, plays essential roles in lipid digestion, absorption, and microbial ecological regulation [27]. Previous studies have also reported inhibitory effects of certain bile acids on Firmicutes taxa [28]. Accordingly, the observed negative correlations are consistent with the regulatory role of bile acids in modulating gut microbial composition. These findings suggest a potential bidirectional interaction in which bile acid metabolites and specific microbial taxa may jointly contribute to the maintenance of intestinal ecological balance.

Overall, this study revealed a multi-level ecological response of the Tibetan sheep intestine to NAPs supplementation. The small intestine exhibited greater sensitivity to exogenous intervention, as reflected by the clear separation of microbial community structure and dynamic modulation of associated metabolic pathways. At the metabolic level, differential enrichment of bile acid, fatty acid, and phospholipid metabolic pathways suggests that NAPs supplementation may influence host energy balance and membrane-associated functions through coordinated regulation of lipid metabolism and bile acid homeostasis [29].Integrated microbiome–metabolome analysis further demonstrated that within the shared variation dimension of bile acid metabolism and lipid-related metabolites [30], specific genera (e.g., Lacrimispora) were significantly associated with multiple metabolites. These associations imply the existence of a potential regulatory network in which microbial metabolism and host-derived metabolic products interact bidirectionally to maintain intestinal ecological balance [31].

It should be noted that these conclusions are derived from statistical associations and inferred biological relevance, therefore carry inherent uncertainty. Functional prediction based on 16S rRNA gene sequencing is constrained by annotation limitations and indirect inference. Although differential metabolite identification and pathway enrichment analyses provide valuable mechanistic clues, causal relationships remain to be validated through targeted metabolite quantification, functional assays, and integration with metagenomic or transcriptomic approaches [32]. Future studies should focus on: (i) high-resolution functional annotation and metabolic pathway reconstruction based on metagenomic or transcriptomic sequencing; (ii) time-series experimental designs to evaluate the dynamic effects and regulatory boundaries of NAPs supplementation; and (iii) mechanistic validation of key microbe–metabolite interactions identified in this study. In summary, this work provides multi-level evidence for understanding the interplay between intestinal microecology and host metabolism in Tibetan sheep and highlights the potential modulatory effects of NAPs supplementation on microbial community structure and metabolic networks across different gastrointestinal sites.

5. Conclusions

This study demonstrates that the small intestine is the primary site of microecological response to NAP supplementation in Tibetan sheep. While the core microbiota remained relatively stable, NAPs significantly reshaped the microbial community structure and diversity (P < 0.05). Integrated omics analysis revealed that these microbial shifts were intrinsically linked to systematic alterations in bile acid, fatty acid, and phospholipid metabolic pathways. Specifically, the strong correlations between key genera and metabolites, such as bile acids and prostaglandins, suggest that NAPs maintain intestinal homeostasis by regulating intestinal microecology and host metabolism. These findings provide a theoretical basis for the application of NAPs as precise, site-specific feed additives to enhance gastrointestinal health and metabolic efficiency in ruminants.