All LC-MS hypergrade solvents (methanol, water, acetonitrile), HPLC grade solvent of methanol, and formic acid were obtained from Merck (Darmstadt, Germany). The standard solution of positive and negative calibration solution (Pierce LTQ positive and negative) for mass spectrometry were purchased from Thermo Fisher Scientific (Rockford, IL, USA).

Beef (n = 5) sample was obtained from five different markets in Yogyakarta and Central Java, Indonesia. Pork (n = 3) and dog meat (n = 3) were purchased from three different meat sellers in Yogyakarta, Indonesia. The loin part was used for metabolomics analysis. The meat samples were stored at −20°C prior used for extraction. Each type of meat was ground using a meat grinder to obtain ground meat, then stored to −20°C until used for analysis. Each meat obtained from all locations were mixed. The samples of pure beef, pork, and dog meat as well as the binary mixtures of beef-pork and beef-dog meat were prepared in triplicates.

The extraction of metabolites for metabolomics analysis was carried out based on the method by Wang et al [14] subjected to slight modifications. Pure BMr samples and those adulterated with DMr and PMr using various concentration levels (0.1%, 1%, 5%, 10%, 25%, and 50% w/w) were prepared. The ground meat was weighed, and the mixtures were prepared in a total weight of 5 g. Samples were placed into a Beaker glass and added with methanol (25 mL). Subsequently, samples were vortexed for 60 s, then ultrasonicated for 30 min at room temperature. The protein was precipitated by storing the samples at −20°C for 1 hour. The supernatant was collected by centrifugation for 10 min at 4°C at 5,000 rpm. The 1 mL of supernatant was pipetted and filtered using PTFE filter 0.22 μm. The filtered supernatant was placed into a clear 2 mL HPLC vial. Three replicates of each sample were prepared.

Metabolomics was performed according to our previous study with slight modifications [25]. Separation of compounds was carried out using an ultra-high performance liquid chromatography (UHPLC Vanquish; Thermo Scientific, Rockford, IL, USA) equipped with a binary pump. Two mobile phases consisted of water containing 0.1% formic acid (A) and methanol containing 0.1% formic acid (B) were used for eluting sample with the flow rate of 0.30 mL/min. An Accucore C-18 column (100 mm×2.1 mm×2.6 μm) was used for separating metabolites in samples which was injected at 10 μL. The temperature of column was set at 40°C during the compound’s separation. Elution was carried out using gradient mode as follows; at first, the mobile phase was set at 95% A, then continue to for 10% A at 16 min. Then, maintained at 10% A until 30 min before setting back to 95% A (35 min). The detection of compounds was performed using a high-resolution mass spectrometer (Orbitrap Q-Exactive; Thermo Scientific, USA). The flow rate setting for sheath gas was 32 arbitrary unit (AU), followed by auxiliary gas of 8 AU and sweep gas of 4 AU, respectively. The setting for parameter of capillary temperature was 320°C followed by gas heater temperature of 30°C. The ionization was performed in both electro spray ionization (ESI) positive and ESI negative modes with MS1 resolution of 70,000 and MS2 resolution of 17,500. The screening of compounds was carried out at m/z range of 66.7 to 1,000 m/z.

The Compound Discoverer (Thermo Scientific, USA) software was used to identify metabolites composition in BMr, DMr, PMr, and their mixtures. The raw total ion chromatogram (TIC) both of sample and blank was used for analysis. Spectrum selection, retention time alignment, feature detection, and metabolite identification were performed. The Chemspider and MzCloud database were used for database matching. Analysis was performed using an untargeted workflow. Only metabolites having full match MS2 fragmentation and error mass between −5 and 5 ppm were selected. The metabolite dataset was converted into an excel file for chemometrics analysis.

The selected metabolites from metabolite identification analysis were used as variables for chemometrics. The data were subjected to unit variance scaling prior to chemometrics analysis. Both pattern recognition and multivariate calibration chemometrics were applied. Chemometrics analysis was conducted using a SIMCA 14.1 (Umetrics, Umea, Sweden) software and Metaboanalyst 5.0 platform. Principal component analysis was performed and observed using PCA score plot, PCA loading plot, R², and Q² values. The PLS-DA was further used and evaluated using R²X, R²Y, Q² values, loading score and loading plot. Validation tests such as permutation test, predicted residual error sum of squares (PRESS) value, and receiver operating characteristics (ROC) test were used for PLS-DA. The identification of discriminating metabolites potential as metabolite markers was conducted using variable importance for projection (VIP) test by identifying the variables with VIP more than 1.0. The metabolites were then subjected to analysis of variance (ANOVA) analysis. Metabolites were with p-value <0.05 were considered as significant metabolites in line with the increase of pork and dog meat levels. Multivariate calibration of PLS as well as its orthogonal form (OPLS) was used to build prediction model. The models were evaluated using the PLS plot, R², and the residual plot. Validation of PLS dan OPLS was performed using root mean square error of cross validation (RMSECV) and root mean square error of estimation (RMSEE).

Metabolomics analysis using LC-HRMS employing an untargeted workflow revealed a wide number of metabolites contained in BMr, DMr, and PMr. Figure 1 shows the TIC from BMr, DMr, and PMr exported from XCalibur software. In general, the samples showed a similar TIC pattern, however, deep investigation found some peak differences at the retention time between 14.0 and 24.0 min. Analysis on the TIC against the databases using Compound Discoverer software could retrieve the metabolite compositions on each sample. Figure 2 illustrates the identified metabolites from beef-dog meat and beef-pork series measured using LC-HRMS metabolomics. Fatty acyls class placed at the highest proportion followed by lipids in both sample series which consist of various types of fatty acids and polar lipids such as phospholipids. Various types of amino acids and organic acids were also identified. The results was in accordance with previous reports as fatty acids, lipids, and amino acids were reported to be present in high amount in meats [26]. In addition, it was found that creatine had the highest peak area in all types of meat (BMr, DMr, and PMr), followed by carnitine. It was in agreement with previous research reporting that meats such as beef, pork, and chicken contained high amount of creatine and carnitine [27]. Due to the complexity of the identified metabolites in BMr, DMr, and PMr samples, therefore, an advanced statistical tool is required to analyze the data for the authentication purposes to guarantee the results.

In this study, PCA was first applied to observe the pattern of samples grouping based on their metabolite compositions. Using seven principal components (PCs), PCA could differentiate pure BMr samples with those adulterated with DMr at several concentration levels, especially in high levels of DMr (10%, 25%, and 50%) with R² = 0.877 and Q² = 0.577. However, in the presence of low concentration of DMr (0.1%, 1%, and 5%), PCA could not clearly differentiate between pure and adulterated samples (Figure 3A). On the other hand, PCA using five PCs subjected to metabolites data from beef-pork series cloud clearly differentiate authentic BMr, BMr adulterated PMr, and pure PMr with R² = 0.680 and Q² = 0.403 as depicted in Figure 3B. The tight cluster of the quality control (QC) samples in the PCA score plot confirmed the stability and the reproducibility of the LC-HRMS method in both sample series. PCA is an unsupervised pattern recognition that has been widely used in many fields of research including metabolomics to find out the similarity and dissimilarity among samples by reducing the number of original variables into few PCs. Samples with score plots appeared close to each other indicates high similarity and vice versa [19]. In this study, the obtained PCA showed that the higher the presence of adulterant in BMr, the closer the score plot to the score plot of adulterant both in beef-dog and beef-pork model. However, the separation among levels in adulterated samples is not clear enough such as a few of score plots are still overlapping. It might be caused by losing some information or data distortion when reducing the data dimensionality [16].

To confirm the PCA results, supervised pattern recognition of PLS-DA was further used to discriminate and classify samples. PLS-DA using four latent variables (LV) clearly discriminated BMr, BMr adulterated DMr, and DMr samples with R²X = 0.797, R²Y = 0.961, and Q² = 0.803. The score plot of PLS-DA in Figure 4A shows that the lowest level of DMr (0.1%) was clearly discriminated from the authentic BMr indicating a good performance of the PLS-DA model for discriminating adulterated samples at low concentration of adulterant. In addition, PLS-DA was also successfully used to discriminate between authentic and adulterated BMr with PMr, with high accuracy as depicted in Figure 4B. Using 9 LV, PLS-DA could detect and discriminate all adulterated samples from the authentic BMr samples with R²X = 0.774, R²Y = 0.960, and R²Y = 0.563. According to the PLS-DA score plot, the higher the level of the adulterants given, the closer the score plot to the adulterant meat, observed in both models. The table of misclassification test showed there is no misclassification found in both PLS-DA model of beef-dog and beef-pork series indicating the high accuracy of the classification models. In addition, the stability and reproducibility of the LC-HRMS method was confirmed by the tight cluster of QC samples.

Supervised PLS-DA is a combination of PLS regression and linear discriminant analysis. The PLS model was first built using the X matrix (independent variables), followed by Y matrix (dependent variables). Subsequently, the LV, also called factors were searched and used for classification. These LV made PLS-DA capable of returning better sample differentiation than using PCA. This result was in agreement with previous studies reporting the suitability of PLS-DA for sample classification of meat from different origins, for detection of pork in halal meats, and for metabolite differentiation of chicken during storage [14,28].

However, supervised pattern recognition such as PLS-DA is susceptible for model overfitting thereby could provide bias result. Therefore, such validation tests including permutation test and ROC test are often required to demonstrate the model validity. The permutation test performed using 999 permutations confirmed the validity of PLS-DA model both in BMr adulterated DMr and BMr adulterated PMr series. In the BMr adulterated DMr model, permutation using six components showed that the original variables had the highest value among all the permutated models with the intersection of Q² value was between zero and below zero, (0.0, −0.231). Meanwhile, in the BMr adulterated PMr, permutation test using 8 components demonstrated the intersection value of Q² of (0.0, −0.616) indicating model validity. The result of ROC test demonstrated the PLS-DA model validity supporting the permutation test. All samples were accurately classified in accordance with the result of misclassification test demonstrated by its AUC (area under the curve) of each class. The AUC value for each sample class obtained from both PLS-DA models was 1, exhibiting a good classification without misclassification [29].

The identification of variables having important role in samples differentiation which associated to discriminating metabolites samples differentiation was performed using the VIP value. These discriminating metabolites were important as potential biomarker candidates to discriminate non-halal meats in beef. The variables having an important role in samples differentiation are shown by the VIP value more than 1.0. The discriminating metabolites from VIP analysis of beef-dog series are shown in Supplementary Table S1. Among those metabolites, 19 metabolites were found to increase with the increase of DMr levels in BMr (Table 1). Meanwhile, the discriminating metabolites from beef-pork series are listed in Supplementary Table S2. Further analysis showed that 28 metabolites were observed to increase with the increase of pork level in beef (Table 2). The discriminating metabolites were cross confirmed using the ANOVA analysis [24]. Metabolites with p<0.05 were significant metabolites, that increase with increasing dog meat and pork levels in beef, respectively. The discriminating metabolites that increased in line with the level of adulterants consisted of various metabolites classes such as amino acids, organic acids, fatty acids, and other lipids. These significant metabolites could be used as the keys to detect and monitor the presence of dog meat and pork, respectively.

On the other hand, the area of metabolites acetyl-L-carnitine, DL-carnitine, C12-carnitine, and C14-carnitine were found to be high in BMr and increased as the level of BMr increased in both series. Carnitine is derived from amino acids, and it is endogenously synthesized in the liver, kidney, and brain from the amino acids of methionine and lysine. Carnitine is reported to be high in beef and it plays an important role in turning fat into energy. In addition, it eliminates some toxic compounds by transporting them out of mitochondria [30]. From the above results, it can be summarized that the utilization of VIP analysis is very beneficial to analyze metabolomics data from untargeted LC-HRMS to obtain discriminating metabolites to discriminate beef from dog and pork meat adulteration.

The metabolites obtained from VIP analysis (VIP>1.0) were subjected to multivariate regression analysis using PLS regression and OPLS regression. Partial least square was successfully used to accurately predict the levels of DMr and PMr in the adulterated beef (% w/w), respectively. A high correlation between actual values and predicted values of DMr in adulterated BMr was obtained (Figure 5A and 5C). The PLS model had a high R² value (0.9986) which is associated to its high model accuracy. In addition, the model had a low value of RMSEE (1.30%) and RMSECV (2.02%), associating to low error and high precision of the PLS model. In addition, PLS also successfully predict the levels of PMr present in BMr with high accuracy and high precision. The model provided R² = 0.9981 with RMSEE of 1.60% and RMSECV of 2.70%. On the other hand, the regression of orthogonal PLS (OPLS) has also been widely applied for predicting target of analytes employing the orthogonal components. In this study, OPLS model was successfully detect and predict the adulterants of DMr and PMr in BMr as exemplified in Figure 5B and 5D. The obtained R² was 0.9986 and 0.9981 with RMSEE of 1.30% and 1.60% for beef-dog and beef-pork series, respectively. In addition, the RMSECV value was 1.92% and 2.73% for beef-dog meat and beef-pork series, respectively.

Analysis on the residual plot both in PLS and OPLS models showed that all the residuals value were randomly spread and located at the standard deviation values between −4 and +4. This result indicated that the residuals were normally distributed, and no outliers were detected. The results showed that both PLS and OPLS regression had a similar performance in predicting targets. The results of PLS and OPLS demonstrated that the discriminating metabolites obtained from VIP analysis were good predictors to build a prediction model to determine levels of dog and pork meat present as adulterants in beef meat. Therefore, PLS and OPLS regression could be used to support the results of VIP analysis in investigating the discriminating metabolites potential as biomarker candidates to differentiate non-halal meats adulteration in halal meats.

PLS-based regression become the most multivariate calibration techniques applied in various analysis to predict the concentration of targeted analytes in the presence of other substances or matrices, because it could be used for accurate prediction of the target. The use of LV to correlate X-matrix and Y-matrix in the regression analysis results in high predicting capacity. Some studies have successfully applied PLS for halal authentication by detection of lard in edible fats and detection of pork in meat and meat-based products using spectroscopic techniques [31]. The discriminating metabolites of pork have been successfully confirmed as good predictors for pork adulteration by PLS and OPLS [25]. Another research also successfully used PLS to verify the role of discriminating metabolites in tea with different storage times. Result reveled that the discriminating metabolites were good predictors of storage time [32]. In our study, the result showed that both PLS and OPLS was potential for accurate prediction of DMr and PMr as the adulterants in halal beef meat.

In metabolomics study, it should be noted that the complexity of metabolomics data and its interpretation is not an easy task. Several factors such as the difference in diets, origins, environmental factors, meat preparation, and storage could affect the metabolites. In this study, we found a panel of metabolites which strongly associated with the increasing levels of dog meat and pork adulteration in minced beef. However, further investigation on the effect of factors that might affect the metabolites is required for future research to verify the metabolites found.