RESULTS

Multivariate analysis of metabolomic profiles

Multivariate statistical analyses were performed to evaluate metabolic differences between HP and LP Sapera goats. PCA was first conducted as an unsupervised approach to visualize the overall variation in serum metabolomic profiles. As shown in Figure 1A, the PCA score plot demonstrated clear separation between the HP and LP groups along the first principal component. The first principal component explained 24.2% of the total variance, whereas the second explained 12.5%. This clustering pattern indicated that the metabolic profiles of the two production groups were distinct.

To further investigate group discrimination, PLS-DA was performed. The PLS-DA score plot (Figure 1B) showed clear separation between HP and LP samples along the first component, indicating distinct metabolic signatures associated with milk production levels. The loading vectors in the biplot indicated that several metabolites contributed strongly to the separation between the two groups.

The predictive performance of the PLS-DA model was evaluated using cross-validation. As illustrated in Figure 1C, the model achieved a high classification accuracy of nearly 1.0. The model fitness reached 1.0, indicating excellent explanation of data variance, whereas the predictive capability reached 0.79, indicating good predictive reliability.

Differential metabolites between HP and LP goats

To identify metabolites associated with differences in milk production, a volcano plot analysis was performed using fold change and statistical significance thresholds. The volcano plot (Figure 2) illustrates the distribution of metabolite features based on log2 fold change and −log10(p-value). Metabolites located on the right side of the plot represent compounds upregulated in the HP group, whereas metabolites on the left side represent compounds enriched in the LP group. Differential metabolites were initially screened for statistical significance and fold change patterns, and the final candidate discriminant metabolites were selected from those consistently contributing to group separation in the multivariate analyses.

Several metabolites showed significant differences between the two production groups (Table 1). Among the candidate discriminant metabolites, N-acetylneuraminic acid (C01352), 2-oxoadipic acid (C04438), fumarate (C00417), pyruvate (C00022), β-hydroxybutyric acid (C05983), glycerol 3-phosphate (C03981), 3-hydroxy-3-methylglutaric acid (C02985), N-carbamoyl-L-aspartate (C00233), and L-aminoadipic acid (C05267) showed higher abundance in the HP group. In contrast, sphinganine (C16566), cystathionine (C05273), phytosphingosine (C15767), sphingosine (C21963), cystathionine-related metabolites (C05269), and ceramide-related intermediates (C21937) were relatively enriched in the LP group.

Notably, sphinganine (C16566) showed the strongest statistical significance (−log10[p] = 9.06), followed by cystathionine (C05273) and N-acetylneuraminic acid (C01352), indicating that these metabolites contributed substantially to the metabolic differences observed between HP and LP Sapera goats. The identified differential metabolites were subsequently used to evaluate biomarkers and metabolic pathways.

Figure 1: Multivariate analysis of metabolomic profiles in HP and LP Sapera goats. (A) PCA score plot showing sample clustering. (B) PLS-DA biplot illustrating group separation and contributing metabolites. (C) Cross-validation performance of the PLS-DA model across components.

Figure 2: Volcano plot of differential metabolites between high-production and low-production Sapera goats.

Key metabolites associated with milk production differences

To identify metabolites that contributed most strongly to discrimination between HP and LP goats, variable importance in projection (VIP) analysis was performed on the PLS-DA model. Metabolites with VIP scores >1.0 were considered influential variables in group separation. The final panel of 15 candidate discriminant metabolites was selected based on their contribution to group separation, statistical significance, and biological relevance to lactation-related metabolic pathways. The top metabolites, ranked by VIP scores, are shown in Figure 3A.

Among the identified metabolites, sphinganine, ceramide-related intermediates, cystathionine, N-acetylneuraminic acid, and fumarate showed the highest VIP scores, indicating that these compounds contributed substantially to the metabolic differences between HP and LP groups. The distribution pattern indicated that several metabolites exhibited contrasting abundance profiles between the two production groups.

To further visualize the relative abundance patterns of these key metabolites across samples, hierarchical clustering analysis was performed and presented as a heatmap (Figure 3B). The heatmap revealed a clear clustering pattern separating HP and LP samples into distinct groups. Metabolites such as N-acetylneuraminic acid, fumarate, β-hydroxybutyric acid, N-carbamoyl-L-aspartate, and 3-hydroxy-3-methylglutaric acid showed higher abundance in the HP group, whereas cystathionine, sphinganine, cystathionine-related metabolites, phytosphingosine, and sphingosine were relatively enriched in the LP group. This clustering pattern confirms the metabolic distinction between HP and LP goats and highlights metabolites that may be associated with differences in milk production performance.

KEGG ID | FC | log2(FC) | raw.pval | −log10(p)
C16566 | 0.3759 | −1.4116 | 8.77E−10 | 9.0569
C05273 | 0.42701 | −1.2276 | 3.26E−07 | 6.4866
C01352 | 2.2732 | 1.1847 | 1.68E−06 | 5.7745
C04438 | 6.7157 | 2.7475 | 0.002927 | 2.5336
C01159 | 0.16053 | −2.6391 | 0.0030748 | 2.5122
C00497 | 0.15012 | −2.7358 | 0.0030942 | 2.5095
C00402 | 3.1749 | 1.6667 | 0.006687 | 2.1748
C05444 | 0.37024 | −1.4335 | 0.009914 | 2.0038
C02220 | 0.34301 | −1.5437 | 0.01294 | 1.8881
C00074 | 0.36457 | −1.4557 | 0.013639 | 1.8652
C05840 | 2.4737 | 1.3067 | 0.013974 | 1.8547
C00033 | 2.5162 | 1.3312 | 0.014119 | 1.8502
C06427 | 2.8646 | 1.5183 | 0.01414 | 1.8496
C01384 | 2.4989 | 1.3213 | 0.014833 | 1.8288
C01056 | 4.3253 | 2.1128 | 0.015631 | 1.8060
C01879 | 2.5763 | 1.3653 | 0.019991 | 1.6992
C00712 | 4.1441 | 2.0510 | 0.022034 | 1.6569
C00134 | 2.4129 | 1.2707 | 0.024447 | 1.6118
C17646 | 2.4057 | 1.2664 | 0.025488 | 1.5937
C01157 | 0.4375 | −1.1927 | 0.02582 | 1.5880
C04487 | 0.4389 | −1.1880 | 0.026551 | 1.5759
C04442 | 0.44449 | −1.1698 | 0.034058 | 1.4678
C12455 | 3.6514 | 1.8685 | 0.035024 | 1.4556
C00089 | 2.2600 | 1.1763 | 0.036372 | 1.4392
C00222 | 3.7416 | 1.9037 | 0.037734 | 1.4233
C01832 | 3.2305 | 1.6918 | 0.042048 | 1.3763
C00455 | 0.40952 | −1.2880 | 0.042816 | 1.3684
C00064 | 3.0001 | 1.5850 | 0.046549 | 1.3321
C00636 | 0.41452 | −1.2705 | 0.047242 | 1.3257
C16533 | 0.26639 | −1.9084 | 0.049568 | 1.3048
C00037 | 4.0666 | 2.0238 | 0.050649 | 1.2954
C00412 | 2.8520 | 1.5120 | 0.053273 | 1.2735
C01013 | 3.6050 | 1.8500 | 0.053288 | 1.2734
C03046 | 2.2774 | 1.1874 | 0.053533 | 1.2714
C05266 | 3.3799 | 1.7570 | 0.069573 | 1.1576
C05161 | 0.26173 | −1.9339 | 0.070993 | 1.1488
C00249 | 2.9188 | 1.5454 | 0.072359 | 1.1405
C16513 | 2.6638 | 1.4135 | 0.074276 | 1.1292
C08323 | 5.2095 | 2.3812 | 0.076216 | 1.1180
C02262 | 3.2420 | 1.6969 | 0.076661 | 1.1154
C05265 | 2.7452 | 1.4569 | 0.078212 | 1.1067

Figure 3: Key metabolites associated with milk production differences. (A) Variable importance in projection scores derived from the PLS-DA model showing the top metabolites contributing to discrimination between HP and LP Sapera goats. (B) Hierarchical clustering heatmap showing the relative abundance patterns of the top metabolites across samples. Red indicates higher relative abundance, whereas blue indicates lower abundance.

Biomarker candidate performance

The relative abundance patterns of selected differential metabolites were further examined to characterize metabolic differences between HP and LP Sapera goats (Figures 4 and 5). Several metabolites, including ceramide-related intermediates, N-acetylneuraminic acid, fumarate, β-hydroxybutyric acid, N-carbamoyl-L-aspartate, pyruvate, 3-hydroxy-3-methylglutaric acid, glycerol 3-phosphate, L-aminoadipic acid, and 2-oxoadipic acid, exhibited higher abundance in the HP group. These metabolites consistently showed elevated normalized intensities in HP samples, indicating their association with metabolic processes related to increased milk production.

In contrast, sphinganine, cystathionine, sphingosine, phytosphingosine, and cystathionine-related metabolites displayed higher abundance in the LP group. The abundance distributions of these metabolites clearly separated LP from HP animals, indicating metabolic alterations associated with lower milk production performance. These differences confirm the presence of distinct metabolic signatures between HP and LP Sapera goats.

Figure 4: Boxplots of differential metabolites with higher abundance in high-production Sapera goats than in low-production Sapera goats.

Figure 5: Boxplots of differential metabolites with higher relative abundance in low-production Sapera goats than in high-production Sapera goats.

Metabolite set enrichment analysis

Metabolite set enrichment analysis was conducted to identify metabolic functions associated with the differential metabolites detected between HP and LP Sapera goats (Figure 6; Table 2). The enrichment results revealed several metabolite sets associated with lipid, energy, and amino acid metabolism.

Figure 6: Metabolite set enrichment analysis of differential metabolites between high-production and low-production Sapera goats

Metabolite set | Total | Expected | Hits | Raw p | Holm p | False discovery rate
Fatty acid elongation | 39 | 0.250 | 3 | 0.00154 | 0.124 | 0.0672
Fatty acid degradation | 39 | 0.257 | 3 | 0.00166 | 0.133 | 0.0672
Citrate cycle (TCA cycle) | 20 | 0.132 | 2 | 0.00696 | 0.550 | 0.188
Glyoxylate and dicarboxylate metabolism | 32 | 0.211 | 2 | 0.0174 | 1.000 | 0.353
Valine, leucine and isoleucine biosynthesis | 8 | 0.0527 | 1 | 0.0516 | 1.000 | 0.836
Fructose and mannose metabolism | 21 | 0.138 | 1 | 0.130 | 1.000 | 1.000
Pyruvate metabolism | 22 | 0.145 | 1 | 0.136 | 1.000 | 1.000
Propanoate metabolism | 22 | 0.145 | 1 | 0.136 | 1.000 | 1.000
Glycolysis / Gluconeogenesis | 24 | 0.158 | 1 | 0.148 | 1.000 | 1.000
Alanine, aspartate and glutamate metabolism | 28 | 0.184 | 1 | 0.170 | 1.000 | 1.000
Lipoic acid metabolism | 28 | 0.184 | 1 | 0.170 | 1.000 | 1.000
Amino sugar and nucleotide sugar metabolism | 31 | 0.204 | 1 | 0.187 | 1.000 | 1.000
Glycine, serine and threonine metabolism | 32 | 0.211 | 1 | 0.192 | 1.000 | 1.000
Cysteine and methionine metabolism | 33 | 0.217 | 1 | 0.198 | 1.000 | 1.000
Arginine and proline metabolism | 35 | 0.230 | 1 | 0.208 | 1.000 | 1.000
Glycerophospholipid metabolism | 36 | 0.237 | 1 | 0.214 | 1.000 | 1.000
Biosynthesis of unsaturated fatty acids | 36 | 0.237 | 1 | 0.214 | 1.000 | 1.000
Valine, leucine and isoleucine degradation | 40 | 0.263 | 1 | 0.235 | 1.000 | 1.000
Tyrosine metabolism | 42 | 0.276 | 1 | 0.245 | 1.000 | 1.000

Among the enriched metabolite sets, fatty acid elongation and fatty acid degradation showed the strongest enrichment signals, each containing three matched metabolites with the lowest raw p-values (p = 0.00154 and p = 0.00166, respectively). These metabolite sets exhibited the highest enrichment ratios, indicating a prominent role of fatty acid metabolism in the metabolic differences between HP and LP animals.

Other enriched metabolite sets included the TCA cycle and glyoxylate and dicarboxylate metabolism, suggesting alterations in central carbon metabolism. Several amino acid-related metabolite sets were also detected, including valine, leucine, and isoleucine biosynthesis; glycine, serine, and threonine metabolism; cysteine and methionine metabolism; and arginine and proline metabolism. Overall, the enrichment analysis indicates that lipid, central energy, and amino acid metabolism contribute to the metabolic differences observed between HP and LP Sapera goats.

KEGG pathway analysis

KEGG pathway analysis integrated pathway enrichment results with pathway topology information, allowing identification of pathways with both statistical significance and biological impact (Figure 7; Table 3). The analysis identified several pathways associated with metabolic differences between HP and LP Sapera goats.

Figure 7: Kyoto Encyclopedia of Genes and Genomes pathway analysis of differential metabolites associated with milk production differences between high-production and low-production Sapera goats.

Pathway name | Match status | p -value | −log(p) | Holm p | False discovery rate | Impact
Fatty acid degradation | 3/39 | 0.0015454 | 2.8110 | 0.12363 | 0.061814 | 0.06358
Fatty acid elongation | 3/39 | 0.0015454 | 2.8110 | 0.12363 | 0.061814 | 0.10791
Citrate cycle (TCA cycle) | 2/20 | 0.0066348 | 2.1782 | 0.51752 | 0.17693 | 0.09637
Glyoxylate and dicarboxylate metabolism | 2/32 | 0.01662 | 1.7794 | 1.00000 | 0.33240 | 0.03000
Valine, leucine and isoleucine biosynthesis | 1/8 | 0.050351 | 1.2980 | 1.00000 | 0.80562 | 0.00000
Fructose and mannose metabolism | 1/20 | 0.12161 | 0.91504 | 1.00000 | 1.00000 | 0.03313
Propanoate metabolism | 1/22 | 0.13300 | 0.87613 | 1.00000 | 1.00000 | 0.04103
Pyruvate metabolism | 1/23 | 0.13865 | 0.85807 | 1.00000 | 1.00000 | 0.19137
Glycolysis or Gluconeogenesis | 1/26 | 0.15540 | 0.80855 | 1.00000 | 1.00000 | 0.09785
Alanine, aspartate and glutamate metabolism | 1/28 | 0.16640 | 0.77885 | 1.00000 | 1.00000 | 0.00000
Lipoic acid metabolism | 1/28 | 0.16640 | 0.77885 | 1.00000 | 1.00000 | 0.00000
Cysteine and methionine metabolism | 1/33 | 0.19334 | 0.71368 | 1.00000 | 1.00000 | 0.00000
Glycine, serine and threonine metabolism | 1/34 | 0.19863 | 0.70195 | 1.00000 | 1.00000 | 0.00000
Arginine and proline metabolism | 1/36 | 0.20913 | 0.67959 | 1.00000 | 1.00000 | 0.00000
Glycerophospholipid metabolism | 1/36 | 0.20913 | 0.67959 | 1.00000 | 1.00000 | 0.01497
Biosynthesis of unsaturated fatty acids | 1/36 | 0.20913 | 0.67959 | 1.00000 | 1.00000 | 0.04545
Valine, leucine and isoleucine degradation | 1/40 | 0.22974 | 0.63876 | 1.00000 | 1.00000 | 0.01084
Tyrosine metabolism | 1/42 | 0.23986 | 0.62003 | 1.00000 | 1.00000 | 0.00000
Amino sugar and nucleotide sugar metabolism | 1/42 | 0.23986 | 0.62003 | 1.00000 | 1.00000 | 0.01063

Among these pathways, fatty acid degradation and fatty acid elongation exhibited the strongest enrichment signals, with the lowest p-values (p = 0.0015) and relatively high pathway impact values of 0.0636 and 0.1079, respectively. These results highlight the important role of lipid metabolism in the metabolic differentiation between HP and LP animals. The matched metabolites contributing to the enriched pathways were interpreted based on KEGG pathway mapping of the selected differential metabolites. Lipid-related pathways, including fatty acid elongation and degradation, were primarily supported by lipid-associated discriminant metabolites, whereas the TCA cycle and carbohydrate metabolism pathways were supported by central carbon metabolism intermediates such as fumarate and pyruvate.

In addition, pathways related to central energy metabolism, including the TCA cycle and glyoxylate and dicarboxylate metabolism, were identified with notable pathway impact values. Other pathways associated with carbohydrate and amino acid metabolism were also detected, including glycolysis or gluconeogenesis, pyruvate metabolism, propionate metabolism, and several amino acid metabolic pathways.

Overall, KEGG pathway topology analysis indicates that lipid metabolism, central carbon metabolism, and amino acid metabolism represent key metabolic processes underlying differences in milk production between HP and LP Sapera goats.