MATERIALS AND METHODS

Ethical approval

All procedures involving animals and biological sample collection were reviewed and approved by the Ethics Committee of the Faculty of Animal Science, Universitas Hasanuddin, Makassar, Indonesia (Approval No. 015/UN4.12/EC/IV/2026). The study was conducted in accordance with institutional guidelines and national regulations governing the ethical use of animals in research. Blood samples were collected by trained personnel using aseptic techniques to minimize pain, stress, and discomfort to the animals. Throughout the experimental period, the goats were maintained under standard management and feeding conditions, and all efforts were made to ensure animal welfare and reduce unnecessary handling during sample collection and processing. No invasive procedures other than routine jugular venipuncture were performed, and no animals were euthanized for the purpose of this study.

Study period and location

The study was conducted from January to April 2026 at the Laboratory of Animal Nutrition, Faculty of Animal Science, Universitas Hasanuddin, Makassar, Indonesia, and the Corpora Science Research Laboratory, Yogyakarta, Indonesia. All laboratory analyses and metabolomic profiling procedures were performed under standardized conditions to minimize environmental variations that could affect metabolic profiles.

Study design

This study was designed to compare and characterize the metabolomic profiles of Sapera goats classified into two production groups. A total of 20 lactating Sapera goats were allocated into two balanced groups comprising 10 HP goats and 10 LP goats. The classification of HP and LP animals, together with animal metadata including parity, days in milk, age, body weight, body condition score, and general health status, followed the records established in our previous study on milk yield efficiency in Sapera dairy goats under tropical conditions [10]. In the present study, emphasis was placed on serum LC-HRMS-based metabolomic profiling and pathway analysis to identify discriminant metabolites and metabolic pathways associated with the HP and LP phenotypes. The classification of animals was based on their milk production performance during lactation. Goats with higher milk yield were categorized into the HP group, whereas those with lower milk yield were assigned to the LP group. All animals were maintained under similar feeding and management conditions to minimize environmental factors that could influence metabolic profiles.

Milk production recording

Milk production was recorded to classify Sapera goats into different production groups. Daily milk yield was measured using a standardized milking procedure [11]. Milk produced by each goat was collected and measured during each milking session using a graduated measuring container, and the total daily milk yield was calculated by summing the milk obtained during morning and afternoon milking.

To obtain representative production values, milk production records were evaluated over three lactation periods. The average milk yield of each animal was calculated across these lactation periods. Animals with consistently higher average milk yield were categorized into the HP group, whereas animals with lower milk yield were classified into the LP group.

Blood sampling and serum preparation

Blood samples were collected in the morning, within the same sampling window and under similar management conditions, before the main feeding period to minimize variation associated with diurnal rhythms and recent feed intake. Samples were collected individually from each goat and were not pooled. Following centrifugation, serum samples were visually inspected, and samples exhibiting obvious hemolysis were excluded from further metabolomic analysis.

Blood samples were collected aseptically from the jugular vein using sterile disposable 21G needles (Terumo®, Tokyo, Japan) and vacuum blood collection tubes without anticoagulant (Vacutainer®, BD, Franklin Lakes, NJ, USA). Blood samples were allowed to clot at 5°C for approximately 30 min, then centrifuged at 3,000 rpm (1,500–2,000 × g) for 10-15 min in a refrigerated centrifuge (Eppendorf® 5810R, Eppendorf, Hamburg, Germany) to separate serum from the clot. The resulting serum was carefully transferred into sterile 1.5-mL microcentrifuge tubes (Eppendorf) using sterile pipette tips (Axygen®, Corning Inc., Corning, NY, USA) [11].

Metabolomic analysis

Serum metabolomic profiling was performed using LC-HRMS following established metabolomics workflows [11]. Metabolites were extracted using a methanol-based protein precipitation method. Briefly, 100 μL of each serum sample was transferred into 1.5-mL microcentrifuge tubes (Eppendorf) and mixed with prechilled 80% methanol (LC-MS grade; Thermo Fisher Scientific™, Waltham, MA, USA) by vigorous vortexing using a vortex mixer (VWR®, Radnor, PA, USA). Samples were incubated on ice for 5 min and centrifuged at 15,000 × g for 20 min at 4°C using a refrigerated microcentrifuge (Eppendorf® 5427R, Eppendorf, Hamburg, Germany). The resulting supernatant was collected and diluted with LC-MS-grade water (Thermo Fisher Scientific) to obtain a final methanol concentration of 53%. The solution was transferred into a fresh tube and centrifuged again at 15,000 × _g _for 20 min at 4°C to remove residual particulates. The final supernatant was transferred into LC-MS vials (Thermo Scientific™, Thermo Fisher Scientific, Waltham, MA, USA) and subjected to LC-MS/MS analysis.

Chromatographic separation was performed using a Thermo Scientific™ Vanquish™ Horizon ultra-high-performance liquid chromatography system (Thermo Fisher Scientific) equipped with a binary pump and an Accucore™ C18 column (100 mm × 2.1 mm internal diameter, 2.6 μm particle size; Thermo Fisher Scientific) maintained at 40°C. The mobile phases consisted of LC-MS-grade water containing 0.1% formic acid and LC-MS-grade acetonitrile containing 0.1% formic acid (Thermo Fisher Scientific). The mobile phase was delivered at a flow rate of 0.3 mL/min over a total run time of 25 min. The gradient elution program started at 5% organic phase, increased linearly to 90% over 16 min, was maintained at 90% for 4 min, and returned to 5% until 25 min for column re-equilibration. The injection volume was 5 μL.

Mass spectrometric detection was carried out using a Thermo Scientific™ Orbitrap™ Exploris 240 high resolution mass spectrometer (Thermo Fisher Scientific) operating in full MS/data-dependent MS2 acquisition mode with polarity switching in both positive and negative electrospray ionization modes. Data acquired in both ionization modes were processed using the same Compound Discoverer workflow. Features from both modes were annotated based on accurate mass, retention time alignment, MS/MS spectral matching, and database-supported identification. When the same metabolite was detected in both ionization modes, annotations with higher confidence, superior MS/MS spectral matching, and stronger database support were prioritized to avoid duplicate metabolite reporting. Full MS scans were acquired at a resolution of 60,000 full width at half maximum across an m/z range of 67-1000 with a mass tolerance of 5 ppm, whereas MS2 spectra were acquired at a resolution of 30,000 full width at half maximum using stepped normalized collision energies of 30, 50, and 70.

Raw LC-MS data were processed using Thermo Scientific™ Compound Discoverer™ software version 3.5 (Thermo Fisher Scientific). Peak detection, retention time alignment, and compound annotation were performed in accordance with standard metabolomics workflows. Raw data and processed features were inspected for peak quality, retention time consistency, mass accuracy, and MS/MS spectral quality. Metabolite identification was carried out by matching MS/MS spectra against the mzCloud spectral library (LC/MS Autoprocessed; LC/MS Reference). Structural prediction was further supported by using the SIRIUS software integrated with the KEGG and PubChem databases. Additional compound annotation was performed using ChemSpider-integrated resources, including the Human Metabolome Database, KEGG, PubChem, Chemical Entities of Biological Interest, ChEMBL, and FooDB. Mass list searches were also conducted using Human Metabolome Database version 5 mass lists and the Natural Products Atlas (release 2024_09) to improve metabolite annotation confidence. Features were retained for further analysis when they exhibited acceptable peak quality, mass accuracy within 5 ppm, and reliable database-supported annotation.

Pooled quality control samples and technical injection replicates were not included in the present LC-HRMS workflow; therefore, quality control-based coefficient of variation filtering and batch-effect correction were not applied.

Bioinformatic analysis

Bioinformatic analysis followed the general workflow described by Yusuf_ et al._ [12]. The processed feature table generated from Thermo Scientific™ Compound Discoverer™ software version 3.5 was used for downstream statistical and bioinformatic analyses. Before statistical interpretation, metabolite features were filtered based on annotation confidence, peak quality, MS/MS matching scores, mass accuracy (≤5 ppm), and database verification against mzCloud, KEGG, PubChem, and the Human Metabolome Database. Features with poor annotation confidence, inconsistent peak detection, or unreliable spectral matching were excluded from further analysis. Identified compounds were annotated using their corresponding KEGG identifiers to ensure standardized metabolite reporting and compatibility with KEGG-based pathway analysis.

The curated metabolite intensity table was imported into MetaboAnalyst 6.0 for statistical analysis. Data were normalized to reduce technical variation among samples, then log-transformed and scaled before multivariate analysis. Missing values, when present, were handled using the missing-value processing function in MetaboAnalyst, in which low-abundance values were replaced with a small-value imputation approach. The normalized dataset was subsequently used for multivariate and univariate statistical analyses.

Comparative metabolomic analyses were performed between HP and LP goats. Principal component analysis (PCA) was applied to evaluate clustering patterns, intrinsic variability, and potential outliers among samples. Subsequently, partial least squares-discriminant analysis (PLS-DA) was performed to identify metabolites contributing to discrimination between HP and LP groups, and model performance was evaluated using coefficients of determination and predictive capability values.

Differential metabolites between groups were identified using a volcano plot analysis that integrated fold change and statistical significance thresholds. Statistical significance was evaluated using Student's t-test or the Mann-Whitney U test when the assumption of normality was not satisfied. Metabolites with p < 0.05 and variable importance in projection scores >1.0 were considered significant contributors to group separation.

Metabolite annotation confidence was reported according to the Metabolomics Standards Initiative framework. Because authentic chemical standards were not used for confirmation, none of the discriminant metabolites were assigned to Level 1. The 15 key discriminant metabolites were classified as putatively annotated metabolites corresponding to Level 2 based on accurate mass, MS/MS spectral matching, and database-supported annotation using mzCloud, KEGG, PubChem, the Human Metabolome Database, ChemSpider-integrated resources, and SIRIUS-assisted structural prediction.

Hierarchical clustering analysis was conducted and visualized using heatmaps to illustrate relative metabolite abundance patterns across samples and groups. Receiver operating characteristic analysis was performed on selected discriminant metabolites to assess their ability to distinguish HP from LP goats. The results were used to support the biomarker potential of the selected metabolites while recognizing that further validation in larger independent populations is required before diagnostic application.

Finally, functional interpretation of differential metabolites was performed using KEGG pathway enrichment and pathway topology analyses in MetaboAnalyst 6.0 to identify metabolic pathways significantly associated with differences in milk production.