*Article* **Fast Quantification of Honey Adulteration with Laser-Induced Breakdown Spectroscopy and Chemometric Methods**

#### **Jiyu Peng 1, Weiyue Xie 1, Jiandong Jiang 1, Zhangfeng Zhao 1, Fei Zhou 2,4,\* and Fei Liu <sup>3</sup>**


Received: 17 February 2020; Accepted: 11 March 2020; Published: 14 March 2020

**Abstract:** Honey adulteration is a major issue in food production, which may reduce the effective components in honey and have a detrimental effect on human health. Herein, laser-induced breakdown spectroscopy (LIBS) combined with chemometric methods was used to fast quantify the adulterant content. Two common types of adulteration, including mixing acacia honey with high fructose corn syrup (HFCS) and rape honey, were quantified with univariate analysis and partial least squares regression (PLSR). In addition, the variable importance was tested with univariable analysis and feature selection methods (genetic algorithm (GA), variable importance in projection (VIP), selectivity ratio (SR)). The results indicated that emissions from Mg II 279.58, 280.30 nm, Mg I 285.25 nm, Ca II 393.37, 396.89 nm, Ca I 422.70 nm, Na I 589.03, 589.64 nm, and K I 766.57, 769.97 nm had compact relationship with adulterant content. Best models for detecting the adulteration ratio of HFCS 55, HFCS 90, and rape honey were achieved by SR-PLSR, VIP-PLSR, and VIP-PLSR, with root-mean-square error (RMSE) of 8.9%, 8.2%, and 4.8%, respectively. This study provided a fast and simple approach for detecting honey adulteration.

**Keywords:** honey; adulteration; feature variable; partial least square regression; laser-induced breakdown spectroscopy

#### **1. Introduction**

Food adulteration is an illegal activity of food production, which may threaten food quality and safety. On one hand, the nutritional value of food is limited because of the reduction of effective components in food. On the other hand, the adulterants may have a detrimental effect on human health. Several scandals concerning food adulteration have been reported around the world [1–3]. Honey is one of the most commonly adulterated foods because of its economical purpose and wide use. There are two main approaches for honey adulteration. One is to mix pure honey with sugar-based adulterants, and the other is to adulterate high-quality honey with inferior honey. These two cases will be explored in this study.

The adulterant usually has a similar constituent or characteristic with the pure honey, and it is hard to distinguish from the appearance. Several studies concerning honey adulteration detection have been reported. Amiry et al. [4] discriminated adulterated honey (mix pure honey with date syrup and invert sugar syrup) with linear discriminant analysis. Different parameters including color indices, rheological, physical, and chemical parameters were used as variables for discrimination. Physical and chemical parameters achieved the best results, with accuracy above 95%. The results highlighted the use of physical and chemical parameters to detect honey adulteration. In addition, Arroyo-Manzanares et al. [5] used gas chromatography-ion mobility spectrometry to detect sugar cane or corn syrup adulterated honey; seven out of nine commercial honeys were classified as adulterated samples. Traditionally, the chemical features of honey are detected with wet chemical analysis, which is time and labor consuming. Hence, several rapid analytical methods based on electronic and optical techniques were proposed by other researchers, e.g., electronic nose [6], electronic tongue [6], fluorescence spectroscopy [7], visible-near infrared spectroscopy [8,9]. The 'fingerprint information' of honey could be rapidly obtained by these sensors, and the adulterated honey could be distinguished with the help of chemometric methods.

For its part, laser-induced breakdown spectroscopy (LIBS), which allows elemental analysis, may be useful for honey authenticity. The elemental information of honey can be obtained through analyzing the atomic emission spectroscopy from plasma which is induced by a laser. It has the advantages of fast detection, multi-elemental analysis, and environmentally friendly feature [10]. As a novel approach in food, it has been used for regional discrimination [11] and elemental detection [12–14]. Because LIBS spectrum often contains numerous variables, chemometric methods are usually used to figure out the useful information and establish models for food adulteration detection. Recently, LIBS was used to classify the botanical origins of honey, and detect rice syrup adulterated samples [15]. However, the adulterant content in honey should be further quantified. Herein, LIBS combined with partial least squares regression was used as an analytical tool for fast quantification of honey adulterant content.

In this study, acacia honey mixed with high fructose corn syrup (HFCS) and rape honey were analyzed by LIBS. The specific objectives were to: (1) analyze the LIBS spectral features of pure honey and adulterants; (2) determine the feature variables that are related to adulteration; (3) quantify the adulterant content with univariate and multivariate analysis.

#### **2. Materials and Methods**

#### *2.1. Sample Preparation*

Honey including acacia honey (Guanshengyuan Co., Ltd, Shanghai, China) and rape honey (Yaoquan Food Co., Ltd, Yunnan, China) were collected from main producers in China, and two kinds of HFCS with different fructose concentrations (F55 and F90) were purchased from markets. HFCS F55 contains 55% fructose, and HFCS F90 contains 90% fructose. In this case, acacia honey was considered as pure honey, and HFCS (F55 and F90) and rape honey were used as adulterants.

Honey adulteration was prepared by mixing the acacia honey with HFCS F55, HFCS F90, and rape honey. To establish models for quantifying adulterant content, acacia honey was adulterated with HFCS and rape honey at 21 different percentages (0%, 5%, 10%, 15%, 20%, 25%, 30%, 35,% 40%, 45%, 50%, 55%, 60%, 65%, 70%, 75%, 80%, 85%, 90%, 95%, and 100%). In addition, adulterated samples for external prediction were prepared at 13 different adulteration rates, i.e., 0%, 8%, 16%, 24%, 32%, 40%, 48%, 56%, 64%, 72%, 80%, 88%, and 96%. The adulteration rates of 0% and 100% indicated pure acacia honey and pure adulterant, respectively. All sample adulteration was performed in three replications, so there were 63 samples for calibration, and 39 samples for prediction. After mixing, all samples were kept in a water bath at 37 ◦C for 12 h to ensure homogeneity.

#### *2.2. LIBS Measurement*

A laboratory-assembled LIBS device was used for honey adulteration detection. The detailed description of the device was introduced in our previous published article [16]. First, 8 g of sample was added in 12-well plates and placed in a X-Y-Z moving stage. A pulse laser (Vlite 200, Beamtech, Beijing, China) operated at 532 nm was used to ablate the sample with energy of 80 mJ. Then, emission light from induced plasma was transferred into an Echelle spectrograph (ME 5000, Andor, Belfast, UK), and detected by an intensified charge coupled device (ICCD, DH334T-18F-03, Andor, Belfast, UK). To improve the signal-to-background ratio, the delay time, integral time, and relative gain of ICCD camera were set at 2 μs, 10 μs, and 26. Single shot scanning was performed in an ablation region of 10 mm × 10 mm with resolution of 1 mm. Hence, 100 successive spectra were collected for each sample, the spectra were averaged to minimize the sample inhomogeneity. Because of the advantages of LIBS, no sample preparation was needed, and the total detection time for one sample was less than two minutes.

#### *2.3. Data Analysis*

Because the peak in LIBS spectrum corresponds to the emission from a certain element or molecule band, the observed peak intensity was used as the variable for analysis. To establish a model for quantifying adulterant content, PLSR was used. In addition, several feature selection methods based on PLSR were used to determine the key LIBS emissions that related to the adulterant content.

PLSR is a commonly and widely used multivariate method for quantitative analysis. It projects the raw variables into new dimensions with the maximal variation, and regresses the first few new variables (latent variable, LV) with respond value [17]. In this case, the raw variables were peak intensities of main emissions, and the respond value was the adulterant content in honey. Before modeling, the auto scale preprocessing method, which used mean-centering followed by dividing each variable by the standard variation of the variable, was used to correct the scaling of each variable. Ten-folds random cross-validation was used to determine the number of LV, and prevent the overfitting. In addition, the straightforward implementation of a statistically inspired modification of the PLS (SIMPLS) algorithm was used to calculate the PLS model parameters [18].

Three feature selection methods including genetic algorithm (GA), variable importance in projection (VIP), and the selectivity ratio (SR) were used in this case. GA is a subset search algorithm that was inspired by biological evolution theory and natural selection [19]. The subset of relevant variables selected by GA is then fitted with PLSR to evaluate the performance, and determine the feature variables. Different from GA, the variable selection based on VIP and SR is carried out by using a threshold of some parameters from the PLSR model. VIP calculates the accumulation of PLS weights, and SR defines the ratio between explained variance and the unexplained variance in the PLS model. The larger values of VIP and SR, the greater contribution of the variable. For the criteria of variable selection, VIP follows the rule of 'greater than one rule', and SR follows the F-test (95%) criterion [20]. In this case, the variables with VIP value greater than 1 and SR value greater than 1.532 were selected as important variables.

After modeling, some measures should be used to evaluate the performance. In this case, model performance was evaluated with correlation coefficient (*r*) and root-mean-square error (RMSE). The *r* value measures the relationship between predicted adulterant content and actual value, and the RMSE value measures the predictive error. The larger the *r* value and the smaller the RMSE value, the better the model performance. All data analyses were carried out in the MATLAB (v2019b, The MathWorks Inc., Natick, MA, USA).

#### **3. Results and Discussion**

#### *3.1. LIBS Spectral Characteristics*

Before quantification, LIBS spectral characteristics of acacia honey, rape honey, HFCS F55, and HFCS F90 were first analyzed (Figure 1). All the LIBS spectra ranged from 240 to 860 nm. In general, the average LIBS spectra for different samples were similar except some emissions in certain spectral range. It was credited to the similar constituent of honey and HFCS. In general, honey contains 75% saccharides (mainly glucose and fructose), 15% water, amino acids, and minerals, etc. HFCS mainly contains glucose and fructose. According to the concentration of fructose, the HFCS can be divided into three categories: F42 (42% fructose), F55 (55% fructose), and F90 (90% fructose). Hence, the main components ablated by laser in both honey and HFCS were glucose and fructose. As shown in Figure 1,

the emissions from C, H, O, and N were observed in all samples. The molecular band CN that usually appears in an organic sample when analyzed in air atmosphere was also found in this case.

**Figure 1.** Average laser-induced breakdown spectroscopy (LIBS) spectrum of honey (acacia honey and rape honey) and high fructose corn syrup (HFCS 55 and HFCS 90).

Some differences in elemental emissions could be observed between honey and HFCS. It was obvious that emissions from Mg, Ca, and K appeared in the spectra of honey, while it cannot be found in the spectra of HFCS. It indicated that the concentrations of Mg, Ca, and K in honey were significantly higher than those in HFCS. In addition, there was no obvious difference between acacia honey and rape honey, except relatively stronger emission of Na in acacia honey. These elemental differences might be used to differentiate the adulterants. However, it was hard to quantify the adulterant content simply by analyzing spectrum. Hence, some modeling methods were further used to quantify the adulterant content.

#### *3.2. Univariate Analysis*

Univariate analysis was used to explore the relationship between adulterant content and single variable and quantify the adulteration. In this case, the peak intensities of main emissions from samples were used for analysis. Univariate analysis was performed by regressing the peak intensity of each emission with the adulterant content, and *r* and RMSE were used to evaluate the results. The corresponding element for each emission could be identified with the National Institute of Standard and Technology (NIST, Gaithersburg, Maryland, USA) database [21]. Table 1 shows the results of univariate analysis between main emission lines and adulterant content. Forty-three univariate models were established. The variables contained emissions from C, Si, Mg, Ca, Na, K, N, H, O, and CN. Four variables with emissions of 748.47, 794.83, 795.17, and 822.43 nm were marked with unknown, because they could not be identified with the NIST database or references.

In general, the models for quantifying adulterant content of HFCS F90 had the best results with higher *r* and lower RMSE. It indicated that high concentration of fructose in HFCS led to greater spectral difference and contributed to the univariate analysis. In addition, for HFCS F90 and HFCS F55, the emissions from Mg II 279.58, 280.30 nm, Mg I 285.25 nm, Ca II 393.37, 396.89 nm, Ca I 422.70 nm, Na I 589.03, 589.63 nm, and K I 766.57, 769.97 nm had compact relationship with the adulterant content, with *r* > 0.9 and RMSE < 11.0%. For rape honey, models based on emissions from Na I 589.03 and 589.63 nm had good results, with *r* of 0.919 and 0.903, and RMSE of 12.0% and 13.0%. It indicated that emissions from mineral elements played an important role in adulteration quantification. It also verified the LIBS spectral difference between acacia honey and adulterants.


**Table 1.** Results of univariate analysis based on peak intensities of main emissions.

Note: The shade color of the table represents the performance of univariate analysis. The shade color of being green indicates the best compact relationship (*r* = ±1) and the lowest predictive error (RMSE = 0).

#### *3.3. Quantification of Adulterant Content Based on Multivariate Analysis*

Multivariate analysis was further used to quantify the adulterant content. First, all variables in univariate analysis were used as the inputs of PLS models. As seen in Table 2, PLS models based on all variables achieved good results for all three types of adulteration. The *r* values for HFCS F55, HFCS F90, rape honey in the prediction set were 0.962, 0.980, 0.988, and the RSME values were 15.6%, 16.6%, 4.7%, respectively. The latent variables for these three models were 4, 4, 5, which were determined by cross validation. The results of PLS models were better than those of univariate analysis. It also verified the advantages of multivariate analysis. The combination of information from multiple emissions contributed to the adulterant content quantification.


**Table 2.** Multivariate analysis results based on partial least square regression (PLSR) and feature selection methods.

Note: No. of LV: number of latent variables; No. of var.: number of variables; C.V.: cross-validation; *r*: correlation coefficient; RMSE: root-mean-square error; GA: genetic algorithm; VIP: variable importance in projection; SR: selectivity ratio.

In addition, results of PLS models based on feature variables (selected by GA, VIP, and SR) are also shown in Table 2. In general, prediction results after feature selection were similar or better than those based all variables. The irrelevant variables in models might worsen the modeling performance [22,23], which also verified the necessity of feature selection. Only one exception happened for the GA-PLS model in HFCS F55 quantification. The RMSE value in prediction set was 0.320, which is greatly worse than that without feature selection (0.156). It might be credited to the selected variables by the GA method. As shown in Figure 2, lots of irrelevant variables were selected. The GA method might not be suitable for feature selection in the honey adulteration with HFCS F55. With the consideration of variable number and prediction performance, the models marked with bold achieved the best results. The RMSE value for HFCS 55, HFCS F90, and rape honey in the prediction set were 8.9%, 8.2%, and 4.8%, respectively. In addition, similar results were achieved in 10-folds cross-validation, and RMSE value for HFCS 55, HFCS F90, and rape honey were 8.5%, 6.5%, and 4.6%, respectively.

We also compared the variables selected with GA, VIP, and SR methods (Figure 2). Row 1, 5, 9 showed the correlation coefficient between each variable and adulterant content of HFCS F55, HFCS F90, and rape honey, respectively. The values of correlation coefficient were in the range of 0 to 1. Other rows represented the variables selected by GA, VIP, and SR methods. Selected variables were represented in blue, and non-selected variables were in white. As shown in Figure 2, VIP and SR methods chose the variables with a high correlation coefficient, while some variables with a low correlation coefficient were selected by the GA method. It was related to the principal of feature selection methods. For the GA method, the variables were randomly combined and verified by PLSR. The variables were selected based on the results of PLSR modeling. For VIP and SR methods, the contribution of each variable was considered in the selection [20]. The variables selected by the GA method might be easily affected when testing with external samples. In addition, VIP and SR methods had some common variables, while the number of selected variables was different. It might be credited

to the different threshold measure of each method. Hence, VIP and SR methods might be recommended for feature selection in quantification of honey adulterant content.

**Figure 2.** Feature variables selected with genetic algorithm (GA), variable importance in projection (VIP), and selectivity ratio (SR) methods. Row 1, 5, 9 shows the univariate analysis result between each variable and adulterant content of HFCS F55, HFCS F90, and rape honey, respectively. Cells with a gradient of blue color indicated the correlation coefficient. Other rows represented the variables selected by GA, VIP, and SR methods. Selected variables were represented in blue, and non-selected variables were in white.

The scatter plot of the best model for quantifying adulteration ratio of HFCS 55, HFCS 90, and rape honey is shown in Figure 3. Among these three models, the quantification for rape honey achieved the best result, with *r* and RMSE of 0.988 and 4.8% in the prediction set. The samples in calibration and prediction sets distributed closely around the regression lines, and the regression lines almost went through original point. The emissions from Mg II 279.58, 280.30 nm, Mg I 285.25 nm, Ca II 393.37, 396.89 nm, Ca I 422.70 nm, Na I 589.03, 589.64 nm, and K I 766.57, 769.97 nm, which were the feature variables in the rape honey quantification, were also included in the other two models. It indicated that these variables might play an important role in honey adulteration analysis.

**Figure 3.** Scatter plot of actual adulterant content vs. LIBS measured adulterant content. Quantification of adulterant content in the mixture of (**a**) acacia honey and HFCS F55 based on the SR-PLSR model; (**b**) acacia honey and HFCS 90 based on the VIP-PLSR model; (**c**) acacia honey and rape honey based on the VIP-PLSR model.

#### **4. Conclusions**

In this study, LIBS combined with chemometric methods was used to detect honey adulteration. The adulterant content of acacia honey (adulterated with HFCS 55, HFCS 90, and rape honey) was successfully quantified. SR and VIP methods detected effectively the most relevant variables for adulteration determination. The emissions from Mg II 279.58, 280.30 nm, Mg I 285.25 nm, Ca II 393.37, 396.89 nm, Ca I 422.70 nm, Na I 589.03, 589.64 nm, and K I 766.57, 769.97 nm were considered as feature variables and played an important role in modeling. The importance of these variables was also verified in univariate analysis. The SR-PLSR, VIP-PLSR, and VIP-PLSR achieved the best results for detecting an adulteration ratio of HFCS F55, HFCS 90, and rape honey, with RMSE of 8.9%, 8.2%, and 4.8%, respectively. The results indicated the promising possibility of using LIBS and chemometric methods for quantification in honey adulteration. In addition, some research concerning model transfer could be explored, and more types of acacia honey as well as adulterants could be included in modeling in further study, which might be helpful for practical application.

**Author Contributions:** Conceptualization F.Z. and F.L.; Data curation, J.P.; Formal analysis, J.P.; Funding acquisition, J.P. and Z.Z.; Investigation, W.X.; Methodology, J.P.; Project administration, F.Z. and F.L.; Resources, J.J. and Z.Z.; Software, W.X.; Validation, W.X., F.Z. and F.L.; Writing – original draft, J.P.; Writing – review & editing, J.P., W.X., J.J., Z.Z., F.Z. and F.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Key R&D Program of China, grant number 2018YFD0700502, Zhejiang Provincial Key Research and Development Program, grant number 2017C02027, And China Postdoctoral Science Foundation, grant number 2019M652143.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Oilseeds from a Brazilian Semi-Arid Region: Edible Potential Regarding the Mineral Composition**

#### **Ivone M. C. Almeida 1, M. Teresa Oliva-Teles 2, Rita C. Alves 1, Joana Santos 1, Roberta S. Pinho 3, Suzene I. Silva 3, Cristina Delerue-Matos <sup>2</sup> and M. Beatriz P. P. Oliveira 1,\***


Received: 17 January 2020; Accepted: 18 February 2020; Published: 21 February 2020

**Abstract:** Oilseeds from five native plant species with edible potential from the Brazilian Caatinga semi-arid region (*Diplopterys pubipetala*, *Barnebya harleyi*, *Croton adamantinus*, *Hippocratea volubilis*, and *Couroupita guianensis*) were investigated regarding their mineral contents. The minerals, Na, K, Ca, Mg, Fe, Cu, Cr, Al, were analyzed by high-resolution continuum source atomic absorption spectrometry (HR–CS AAS) and P by the vanadomolybdophosphoric acid colorimetric method. K, Mg, and P were the main elements found (1.62–3.7 mg/g, 362–586 μg/g, and 224–499 μg/g dry weight (dw), respectively). *B. harley* seeds contained the highest amounts of K and P, while *C. guianensis* seeds were the richest in Mg. Fe was the most abundant oligoelement (2.3–25.6 μg/g dw). Cr contents were below the limit of quantification for all samples and Al amounts were low: 0.04–1.80 μg/g dw. A linear discriminant analysis clearly differentiated *B. harleyi* and *C. guianensis* samples from the remaining ones. In sum, these oilseeds from the Brazilian Caatinga semi-arid region seem to have the potential to be used as natural sources of minerals, mainly K.

**Keywords:** oilseeds; Caatinga; native; minerals; spectrometry

#### **1. Introduction**

The plant species that grow in the arid and semi-arid land regions of the world have attracted considerable attention in recent years. The Caatinga, the main ecosystem in the Northeastern Brazil (~800.000 km2), is a seasonally dry tropical forest composed by a heterogeneous mix of plant species, primarily deciduous xerophytic spiny shrubs, and trees [1,2]. Despite the huge diversity, only a few of the species of this semi-arid region are exploited for industrial purposes [2]. The rural populations living in this area use many of these plants for food, fuel, timber, medicines, and livestock feed. These products contribute to income generation and offer an alternative source of livelihood when crops fail under the erratic and low rainfall conditions. However, human intervention over the years, such as deforestation, livestock pasture, and cropping, has contributed to the loss of vegetation cover and biodiversity, and extensive soil degradation in this biome [3–9].

Although different studies have already identified some plant species with great economic potential [2,3,8–11], the lack of scientific knowledge for the majority remains. The evaluation of the chemical composition of non-conventional tropical plants is essential not only for their fundamental knowledge, but also to support their use as raw material for food, chemical, and pharmaceutical industries [5]. For instance, in a recent study, Andrade et al. showed that native species from Caatinga used by the local population to treat inflammatory disorders have also a good photoprotective potential and could be used for pharmaceutical preparations [11]. In addition, nuts and seeds from the tropics have been shown to be rich in lipids and relatively good sources of inexpensive and renewable carbohydrates and proteins, and compounds with high added value, like lipid-soluble vitamins, phytosterols, and other phytochemicals, as well as notable quantities of minerals [12–14].

The insufficiency of minerals in humans is mainly caused by an unbalanced diet. In fact, Fe, Ca, Mg, and Cu are among the mineral elements most commonly lacking in human diets [15], and their deficiency could have a negative impact on health. The information about the mineral composition of unconventional oilseeds from Caatinga is very scarce. The aim of this work was, therefore, to analyze the mineral profile of oilseeds from five native species of this region using high-resolution continuum source atomic absorption spectrometry (HR–CS AAS), in order to evaluate their edible potential as a source of minerals. The seeds selected for this study are from Brazilian plants traditionally used in folk medicine, namely, *Croton adamantinus* Mull. Arg. (Euphorbiaceae), *Barnebya harleyi* W.R. Anderson & B. Gates *Diplopterys pubipetala* (A. Juss.) W.R. Anderson & C. Davis, (Malpighiaceae), *Hippocratea volubilis* L. (Celastraceae), and *Couroupita guianensis* Aubl. (Lecythidaceae). The present work is, as far as we know, the first published report on the mineral composition of these oilseeds.

#### **2. Materials and Methods**

#### *2.1. Chemicals and Standard Solutions*

Ultrapure water from a Simplicity 185 system (resistivity 18.2 MΩ cm; Millipore, Belford, USA) and nitric acid (65%; Suprapure®, Merck, Darmstadt, Germany) were used for the preparation of samples and standards. All reagents were of analytical or suprapure reagent-grade. Working standard solutions were prepared from the stock solutions of Na, K, Ca, Mg, Fe, Cu, Cr, and Al (1000 mg/L; Panreac, Barcelona, Spain) by proper dilution with a 0.5% and 1% (*v*/*v*) nitric acid solution. Cesium chloride (0.1% *w*/*v*; Sigma-Aldrich, Steinheim, Germany) was used as flame ionization chemical suppressor. For Ca analyses, lanthanum chloride heptahydrate (1%; *w*/*v*) was added to all solutions in order to minimize the oxy–salts interferences in FAAS analysis. A 0.1% (*w*/*v*) amount of magnesium nitrate (10 g/L; Panreac, Barcelona, Spain) was used as matrix modifier for the electrothermal atomic absorption spectrometry (ETAAS) quantification of Fe and Al; 0.05% (*w*/*v*) magnesium and 0.1% (*w*/*v*) palladium nitrates (10 g/L; Merck, Darmstadt, Germany) for Cu analysis.

Potassium dihydrogen phosphate (99.5%; Riedel-de Haën, Seelze, Germany), ammonium heptamolybdate tetrahydrated (99.0%, Merck, Darmstadt, Germany), and ammonium monovanadate (99.0%, Merck, Darmstadt, Germany) were the reagents used for P analysis. Working standards solutions were manually prepared for HR–CS FAAS and P analysis. Solutions for HR–CS ETAAS calibration curves were automatically prepared by the autosampler.

All glassware and plastic materials were cleaned by treatment with nitric acid (10%, *v*/*v*) for 24 h and then rinsed with water, prior to use.

#### *2.2. Samples*

Mature fruits of five different indigenous plant species (*Diplopterys pubipetala*, *Barnebya harleyi*, *Croton adamantinus*, *Hippocratea volubilis*, and *Couroupita guianensis*) were obtained from five individuals of each species in Pernambuco State (Brazil).

The fruits were dehydrated at 60 ◦C for 48 h, until constant weight, and seeds were manually removed from the fruits. Seeds of each species were thinly ground in a mill and stored in plastic containers at 4 ◦C until use.

#### *2.3. Sample Preparation*

Aliquots (1 g) of the dried ground seed specimens were dry-ashed in a muffle furnace at 500 ◦C, overnight. To the resulting white ash, 200 μL of nitric acid (65%) was added and ashed for another hour at 500 ◦C. An amount of 5.0 mL of nitric acid (65%) was used to re-dissolve minerals and, afterward, evaporated to a final volume of 2 mL. After filtration, the solutions were diluted to 50 mL with ultrapure water. Blanks were also carried out. The mineralization procedure was done in triplicate for each composite sample [16].

#### *2.4. Analytical Methodologies*

An Analytik Jena contrAA 700 atomic absorption spectrometer (Analytik Jena, Germany), equipped with a high-intensity xenon short-arc lamp as continuum radiation source, and automatic MPE60 and AS52s autosamplers (Analytik, Jena, Germany), for electrothermal and flame atomization, respectively, was used throughout this work. Air-acetylene flame (high-purity grade) was used for K, Mg, Ca, and Na quantification at 766.4908 nm, 285.2125 nm, 422.6728 nm, and 588.9953 nm, respectively. Transversely heated graphite furnace with platform-pyrolytically coated graphite tubes (Analytik Jena, Germany) and high-purity argon was used for electrothermal determinations of Fe (248.3270 nm), Cu (324.7540 nm), Cr (357.8687 nm), and Al (309.2713 nm). Pyrolysis and atomization temperatures were optimized to maximize absorbance signals and to minimize backgrounds and matrix interferences. The optimized electrothermal programs used are shown in Table 1. Measurements were performed at the selected primary atomic lines at a wavelength-integrated absorbance equivalent to 3 pixels (central pixel ± 1).

**Table 1.** Temperature programs for Fe, Cu, Cr, and Al determination in seeds by high-resolution continuum source electrothermal atomic absorption spectrometry (HRCS ETAAS).


The method accuracy was assessed through recovery experiments of the seed samples spiked with three different concentrations (50%, 100%, and 125% of the sample contents) of the aqueous standards. Precision was evaluated through intra- and inter-day measurements, by performing three sample measurements on the same day (intra-day) and over three days (inter-day). The limits of quantification (LOQ) were based on the residual standard deviation (10-fold) of calibration curves.

*P* was determined by the vanadomolybdophosphoric acid colorimetric method according to the Standard Method 4500-P [16]. The color development reagent for P determination was prepared by the addition of ammonium heptamolybdate tetrahydrated and ammonium monovanadate. Phosphorus measurements were performed at 420 nm in a dual-beam UV/Vis spectrophotometer (Evolution 300, Thermo Scientific, Boston, MA, USA).

#### *2.5. Statistical Analysis*

The results were reported as mean ± standard deviation. One-way analysis of variance (ANOVA) and Tukey's HSD post-hoc test were applied to evaluate the possible significant differences (*p* < 0.05) in mineral content between the different seeds. The mineral compositions of the seeds were also compared using multivariate statistical analysis. A linear discriminant analysis (LDA) was performed with forward stepwise analysis, considering a *p*-value of 0.01, using Wilk's lambda (λ) as the selection criterion, and an F-statistic value to determine the significance of the changes λ when a new variable

is tested. Followed this was canonical analysis to retrieve more information about the nature of the discrimination between the samples. Statistical analysis was carried out using STATISTICA software v.10 for Windows (Statsoft Inc., Tulsa, USA).

#### **3. Results and Discussion**

In previous research, Pinho et al. already identified *D. pubipetala*, *B. harleyi*, *C. adamantinus*, and *H. volubilis* as promising species for cultivation, based on their oil content and fatty acids profile [13]. *C. adamantinus* is used in the treatment of wounds, due to its anti-inflammatory and analgesic properties. Their seeds contain 37.1% of oil, with high concentrations of essential linoleic (44.2%) and linolenic (45.2%) fatty acids, and are rich in vitamin E, interesting characteristics for the human diet [13,17]. *B. harleyi* and *D. pubipetala* seeds contain high amounts of protein and are promising sources of oil (46.4% and 45.3%, respectively), in particular, of α-linolenic acid (42.8% and 31.9%, respectively) [13]. *H. volubilis* is a traditional expectorant. The seeds contain edible oil (45%–50%), but with a slightly bitter taste. *C. guianensis* leaves, flowers, and barks are known to contain many active ingredients that possess antibiotic, antifungal, antiseptic, antinociceptive, and immunostimulant activities ([18–20]. They contain about ~30% of oil, consisting mainly in linoleic acid (>80%) [13,21,22].

The minerals analyzed in this work were selected based on their physiological importance (e.g., Ca, Mg, K, ... ) or relevance as eventual indicators of environmental contamination (e.g., Cr and Al), since the perception of their composition is essential to better understand their potential to be used in human and/or animal nutrition or, for instance, in the dietary supplement industry.

Flame (FAAS) and electrothermal atomic absorption spectrometry (ETAAS), inductively coupled plasma optical emission spectrometry (ICP–OES), and inductively coupled plasma mass spectrometry (ICP–MS) are the main methods employed in mineral analysis [23]. In this work, the more recent technique HR–CS AAS was employed to evaluate the mineral content of the referred seeds. This equipment presents main advantages over the traditional line source AAS, namely the visualization of the spectral environment in the vicinity of the analytical line, the simultaneous background correction, and the improvement in precision and detection limits due to the higher signal-to-noise ratio [24,25].

#### *3.1. Method Validation*

Method parameters including linear range and correlation coefficient (*r*), LOQ values, repeatability, and mean recovery percentages are presented in Table 2. Na, K, Ca, Mg, P, Fe, Cu, Cr, and Al were quantified in seeds by the external calibration method with aqueous standard solutions. Good linearity was achieved for all elements over the range of tested concentrations, with correlation coefficients higher than 0.997 for all calibration plots. The LOQ ranged from 4.6 (Ca) and 11 μg/g (Ca) for major elements, and from 0.025 (Cu) to 0.11 μg/g (Al) for minor elements. Accuracy and precision of the method were good: recoveries ranged from 85% to 114%, intra-day precisions from 1.5% to 5.2%, and inter-day relative standard deviation (RSD) between 2.1% and 7.9%. These data confirmed the suitability of the method for the designated application.


**Table 2.** Methods validation parameters: linear range and correlation coefficients (*r*) of the analytical calibration curves, intra- and inter-day precisions (RSD, %), and limits of quantification (LOQ) for the elements studied.

<sup>a</sup> Average RSD: inter-day (*n* = 3); intra-day (*n* = 3, three days).; <sup>b</sup> Mean recovery of experiments with three levels of sample spiking (50%, 100%, and 125% of the sample contents).

#### *3.2. Seed Analysis*

The potential of the selected oilseeds as sources of minerals (Na, K, Ca, Mg, P, Fe, Cu, Cr) was ascertained and is described in Table 3.


**Table 3.** Mineral concentrations (μg/g dry weight (dw)) of the analyzed oilseeds.

Values are presented as mean ± standard deviation of triplicate determinations. Within each line, different letters represent significant differences between species, at *p* < 0.05.

Physiologically, the macroelements Na, K, Ca, and Mg are essential for several human body functions, while P is the principal reservoir for metabolic energy and a co-factor of many enzymes [26]. K was the most prevalent element in the mineral composition of these seeds, followed, in decreasing order, by Mg, P, Ca, and Na. Among the studied samples, *B. harleyi* contained the highest K concentration (3649 μg/g). K contents ranged from 1636 to 3649 μg/g, levels considerably higher than those reported by Onyeike and Acheru [27] for castor seeds, coconut seeds, dikanut seeds, groundnut seeds, melon seeds, oil bean seeds, and palm kernel seeds (134–168 μg/g). This is an interesting characteristic taking into account the health importance of K.

The Mg contents in the samples (362 to 586 μg/g) were very modest compared to the levels found for some common culinary nuts and seeds (almonds (~2000 μg/g), pecans (~1500 μg/g), hazelnuts (2000 μg/g) [28], and mustard oilseeds (~7000 μg/g) [29].

The P content ranged from 224 to 499 μg/g. These results are lower than those reported for umbu (*Spondias tuberosa* Arr. Cam) seeds (7720–8250 μg/g) [10] and baru (*Dipteryx alata* Vog.) almonds (7300 μg/g) [30], two Brazilian plant species.

Ca levels (134 to 227 μg/g) were lower than those reported for pumpkin seeds (390–420 μg/g), sunflower seeds, and pistachios (~900 μg/g), and peanuts (520–530 μg/g) [28], but seeds from *D. pubipetala*, *B. harleyi*, and *H. volubilis* presented Ca values identical to pine nuts (130 μg/g) [28]. The richest source of Ca were the seeds of *C. guianensis* (227 μg/g), with levels slightly higher to those found in the nuts of *Cyperus esculentus* (188 μg/g) [31].

The Na content was low for all the samples (<106 μg/g), another good characteristic for edible seeds having in view its health effects.

Fe and Cu are essential trace elements. About two-thirds of the total body iron are involved in metabolic or enzymatic functions, mostly in erythrocytes (hemoglobin) and muscles (myoglobin) while the physiological functions of Cu arise directly from its role in a number of Cu containing enzymes, being essential for brain growth and maturation. Cr is also an essential element usually present in food in the trivalent form. In living organisms, Cr (III) is a component of enzymes that controls glucose metabolism and synthesis of fatty acids and cholesterol. Cr (VI) is toxic and carcinogenic to humans, but it is not usually present in food [32,33].

From a food safety perspective, Al contents of the oilseeds were also evaluated. Most of the population is not at risk for Al toxicity since humans have diverse mechanisms to prevent significant absorption and to aid its elimination. However, when protective gastrointestinal mechanisms are bypassed, renal function is impaired or exposure is high, Al can accumulate in the body, and originate impaired neurological development, Alzheimer's disease, metabolic bone disease, dyslipidemia or even genotoxic activity [34].

The levels obtained for the minor minerals Cu, Cr, and Al were generally low for all samples. Fe was the most abundant of the analyzed oligoelements in all seeds, ranging from 2.34 μg/g in *C. adamantinus* to ten-fold more (25.6 μg/g) in *B. harleyi*. These values are lower than those described for mustard seeds (75–170 μg/g) [29], but the Fe level in *B. harleyi* is comparable to the values found in the literature for some well-known edible nuts such as pecan nuts (24–26 μg/g) and peanuts (20 μg/g) [28], and to the tropical babassu (*Orbignya speciosa*) nuts (18–33 μg/g) and sapucaia (*Lecythis pisonis*) nuts (21–36 μg/g) [14]. *H. volubilis* and *C. guianensis* oilseeds contained similar Fe concentrations to those reported for cupuassu (*Theobroma grandiflorum*) seeds (7.5 μg/g) [14]. *D. pubipetala* and *H. volubilis* seeds contained less than quantifiable levels of Cu (<0.025 μg/g) and the remaining samples contained lower amounts than those found in mustard seeds (6–13 μg/g). Cr concentrations fell below its quantification limit (<0.039 μg/g) in all seeds. Aluminum concentrations ranged from 0.044 μg/g in *D. pubipetala* to 1.80 μg/g in *C. guianensis* seeds, less than the concentration indicated for *Allantoma lineata* seeds (8 μg/g), an Amazon tree from the same family [22].

In sum, *B. harley* seeds had the highest content of K, P, Na, Fe, and Cu, while *C. guianensis* presented the highest levels of Mg, Ca, and Al from the five analyzed seeds.

Although all the oilseeds analyzed showed, in general, a similar mineral profile (high concentration of K followed by Mg > P > Ca > Na > Fe > Al > Cu), their composition showed several significant (*p* < 0.05) differences between them (Table 3). A Linear Discriminant Analysis (LDA) was performed to determine which minerals were able to discriminate between two or more naturally occurring groups (Table 4).


**Table 4.** Summary of discriminant function analysis.

Wilks' Lambda: 0.00000 approx. F (28.15) = 56.326, *p* < 0.0000.

The model was generated through forward stepwise analysis, considering variables one by one, choosing those with the most significant contribution to the data discrimination. Cr content was not considered in this multivariate analysis, as it was always inferior to the LOQ of the analytical method. All the other minerals (K, Mg, P, Ca, Na, Fe, Cu, and Al) were selected for the model generated by the LDA, Al and Fe being the variables that showed the highest contribution to the model. A subsequent canonical analysis generated, in total, four canonical variates (CV). The minerals and the case samples were displayed in a canonical variate scatterplot (Figure 1) of the first two discriminant functions, which comprised 98.9% of the information contained in the generated model. The first dimension (CV 1) represents 70.0% of the data variance, being more related to their Al, Fe, and Cu content. The second canonical variate (CV 2) adds more than 28.9% of information to the model, separating the samples according to their Fe, K, Ca, and Mg content. The relations of the variables with the first two canonical variates are presented by the resulting equations:

CV 1 = −4.18[Al] + 3.34[Fe] − 0.97[Mg] − 0.048[P] + 1.29[Ca] − 2.65[Cu] − 1.73[K] − 0.81[Na] (1)

$$\text{CV} \, 2 = -0.15 \text{[Al]} + 1.19 \text{[Fe]} - 0.75 \text{[Mg]} + 0.13 \text{[P]} + 0.83 \text{[Ca]} + 1.06 \text{[Cu]} + 1.23 \text{[K]} + 0.43 \text{[Na]} \tag{2}$$

**Figure 1.** Scatterplot of all analyzed cases on the first canonical variate (CV1) versus the second canonical variate (CV2) (A: *Diplopterys pubipetala*; B: *Barnebya harleyi*; C: *Croton adamantinus*; D: *Hippocratea volubilis*; E: *Couroupita guianensis*).

From the graphical display of the studied cases (see Figure 1), there were two seed species, namely *B. harleyi* and *C. guianensis*, visibly separated from the others. In the CV1, the Al content of the samples allowed to clearly distinguish the separation between the *C. guianensis* seeds and the others. The *B. harleyi* seeds were discriminated by their highest content of Fe, K, and Al content. The seeds *D. pubipetala*, *C. adamantinus*, and *H. volubilis* were displayed more closely to the positive side of the CV1 showing that the mineral composition of these three seeds was more similar.

#### **4. Conclusions**

The mineral contents of oilseeds from five native plant species from Brazil were evaluated. All the samples analyzed seem to be good sources of K, Mg, and P and presenting also low levels of Na. Nevertheless, *B. harley* oilseeds presented significantly higher amounts (*p* < 0.05) of K and P, Fe, and Cu compared with the other species. *C. guianensis* seeds were the best source of Mg. These two samples were, indeed, clearly differentiated from the remaining ones through linear discriminant analysis. In general, Cr and Al contents were low for all samples. The results show that these oilseeds from the Brazilian Caatinga semi-arid region have the potential to be used as food products and as mineral sources (mainly K).

**Author Contributions:** Conceptualization, R.S.P., S.I.S. and M.B.P.P.O.; Formal analysis, I.M.C.A., R.C.A. and J.S.; Funding acquisition, C.D.-M., M.B.P.P.O.; Investigation, I.M.C.A. and M.T.O.-T.; Methodology, M.T.O.-T.; Project administration, S.I.S. and M.B.P.P.O.Resources, R.S.P., S.I.S. and C.D.-M.; Supervision, M.T.O.-T., C.D.-M. and M.B.P.P.O.; Validation, M.T.O.-T.; Writing-original draft, I.M.C.A., R.C.A. and J.S.; Writing-review & editing, M.T.O.-T. and R.C.A.

**Funding:** The work was supported by UIDB/50006/2020 with funding from FCT/MCTES through national funds.

**Acknowledgments:** R. C. Alves is grateful to *Fundação para a Ciência e a Tecnologia* for the CEECIND/01120/2017 contract.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Rapid Identification and Visualization of Jowl Meat Adulteration in Pork Using Hyperspectral Imaging**

#### **Hongzhe Jiang \*, Fengna Cheng and Minghong Shi**

College of Mechanical and Electronic Engineering, Nanjing Forestry University, Nanjing 210037, China; cfn1218@163.com (F.C.); mhshi@njfu.edu.cn (M.S.)

**\*** Correspondence: jianghongzhe@njfu.edu.cn

Received: 31 December 2019; Accepted: 4 February 2020; Published: 6 February 2020

**Abstract:** Minced pork jowl meat, also called the sticking-piece, is commonly used to be adulterated in minced pork, which influences the overall product quality and safety. In this study, hyperspectral imaging (HSI) methodology was proposed to identify and visualize this kind of meat adulteration. A total of 176 hyperspectral images were acquired from adulterated meat samples in the range of 0%–100% (w/w) at 10% increments using a visible and near-infrared (400–1000 nm) HSI system in reflectance mode. Mean spectra were extracted from the regions of interests (ROIs) and represented each sample accordingly. The performance comparison of established partial least square regression (PLSR) models showed that spectra pretreated by standard normal variate (SNV) performed best with Rp <sup>2</sup> = 0.9549 and residual predictive deviation (RPD) = 4.54. Furthermore, functional wavelengths related to adulteration identification were individually selected using methods of principal component (PC) loadings, two-dimensional correlation spectroscopy (2D-COS), and regression coefficients (RC). After that, the multispectral RC-PLSR model exhibited the most satisfactory results in prediction set that Rp <sup>2</sup> was 0.9063, RPD was 2.30, and the limit of detection (LOD) was 6.50%. Spatial distribution was visualized based on the preferred model, and adulteration levels were clearly discernible. Lastly, the visualization was further verified that prediction results well matched the known distribution in samples. Overall, HSI was tested to be a promising methodology for detecting and visualizing minced jowl meat in pork.

**Keywords:** hyperspectral imaging; jowl meat; minced pork; meat adulteration; visualization

#### **1. Introduction**

Meat always plays an important role in the constitution of human diets around the world. People consume meat mainly due to its rich nutritional contents of essential amino acid and vitamins [1]. With the continuous and rapid increase of meat consumption nowadays, meat quality and safety control has become the major priority [2]. Meat is quite susceptible to suffering adulteration, such as the known horsemeat scandal in 2013 that horsemeat was detected in beef. From the perspective of consumers, they highly desire reliable and clear information about the meat or meat products they purchase [3]. Therefore, regardless of deliberate, accidental, or economically motivated meat adulteration, rapid identification is always one of the main issues in further prevention.

Minced meat, also named ground meat, is one of the most popular meat types. It is versatile in that it can be the major ingredient of a variety of meat products including sausages, patties, hamburgers, and meatballs [4]. Especially in Chinese diet culture, traditional dumplings, steamed stuffed buns, or wontons fillings are also mainly composed of this important ingredient. However, the morphological structure of muscles is removed when meat is minced so that the occurrence of adulteration in minced meat is hardly recognized by visual analysis. Pork jowl meat is a kind of lymphatic meat, which contains thyroid, lipomas, and significant amounts of lymph nodes [5]. It is low-price, full of stench, and should be discarded after the animals are slaughtered [5]. In the present period, a relatively high rate of fraudulent phenomena occur in China involving minced jowl meat being substituted for or added into minced pork to be served as fillings in food. These incidents seriously pose health treats and violate rights for consumers. Meat adulteration is not routinely detected, thus, it is highly desirable to rapidly identify if the fillings are adulterated with jowl meat.

Several detection techniques including high-performance liquid chromatography [6], polymerase chain reaction [7], mass spectrometry [8], differential scanning calorimetry [9], and enzyme-linked immunosorbent assays [10] were demonstrated to be effective in detecting meat adulteration. However, these techniques need complex sample preparation and destructive operation, which are high-cost, time-consuming, and laborious. More recently, a considerable number of fast optical techniques have shown potential in detecting meat adulteration. Among them, UV-visible [11], near-infrared spectroscopy (NIRS) [12,13], mid-infrared spectroscopy [14], Raman spectroscopy [15], laser induced breakdown spectroscopy [16], and computer vision [17] have been successfully developed to identify various adulterants in meat. Typically, conventional NIRS could overcome the above-mentioned drawbacks of traditional techniques for its sensitive, rapid, reagent-free, and nondestructive nature, and it has great potential to adapt real-time monitoring application [18,19]. However, the main limitation is that the spectral single-point detection in preselected areas cannot well represent the whole heterogeneous meat sample. In this light, interest in hyperspectral imaging (HSI) is continuously growing. HSI integrates conventional spectroscopy and imaging to acquire both spectral and spatial information from an identical sample to provide the traceable chemical and physical qualities simultaneously [20]. It has received a lot of attention and interest in inspecting both raw and processed meat items [21,22].

With regard to the detection of meat adulteration using HSI, several examples have previously been reported, for instance, pork or beef adulterated with chicken [17], chicken adulterated with carrageenan [23], beef adulterated with pork [24,25], prawn adulterated with gelatin [26], lamb adulterated with pork [27], beef adulterated with horsemeat [28,29], lamb adulterated with duck [30], etc. All the above studies achieved good performance through a combination of HSI and chemometrics. However, to date, few studies focused on the adulteration of minced pork with minced jowl meat.

Therefore, the main aim of this study is to investigate the feasibility of visible and near-infrared (VNIR, 400–1000 nm) HSI for detecting minced jowl meat in pork. The specific objectives were (1) to establish partial least square regression (PLSR) models based on spectra extracted from hyperspectral images and (2) to identify effective wavelengths for developing multispectral models and visualizing the adulteration.

#### **2. Materials and Methods**

#### *2.1. Sample Preparation*

Pure pork meat (from *Longissimus dorsi* muscle) and adulterant of jowl meat from the homologous slaughtered porcine body were purchased from a local retail market in Nanjing, China and transported to our laboratory within 30 min. On the day of purchase, pork meat was first cut into small pieces and minced using a meat grinder (S2-A808, Joyoung Co. Ltd., Jinan, China) for 60 s. The grinder was carefully washed using detergent and hot water, rinsed with distilled water, wiped with paper towel, and totally dried before jowl meat mincing use. In the same way, jowl meat was subsequently minced to be used as the adulterant. Minced pork samples were adulterated by mixing the minced adulterant into pork in range of 10%–90% (w/w) at 10% increments. Additionally, pure pork (0% adulteration level) and pure jowl meat (100% adulteration level) samples were also prepared. Sixteen replicates were prepared for each adulteration level, and samples were individually weighted. They were thoroughly mixed together to obtain a roughly homogenous paste with a final weight constant at 50 g. Then, all the prepared samples were put into round disposable Petri dishes (9 cm in diameter × 1.4 cm deep) with flat surfaces for subsequent imaging procedure.

In total, 176 samples including 16 samples at each adulteration level (11 levels) were prepared. Among them, three quarters of the samples at each adulteration level (12 samples × 11 levels = 132 samples) were randomly selected to be assigned to calibration set, and the residual one quarter (4 samples × 11 levels = 44 samples) were used for purpose of independent prediction. Furthermore, in order to prove the creditability of visualization results, two more control samples with known distributed patterns were also prepared. As for control sample one, four fan-shaped parts at different individual adulteration levels of 100%, 80%, 40%, and 20% were included. With regard to control sample two, two semicircular areas at individual 0% and 50% adulteration levels were covered. After that, the hyperspectral images of all the prepared samples were subsequently captured.

#### *2.2. Hyperspectral Image Acquisition and Calibration*

All the prepared samples were scanned using a laboratory-based push-broom hyperspectral imaging system in reflectance mode. The system converted the visible and near-infrared (VNIR) spectral range of 400–1000 nm (284 spectral bands) to capture the hyperspectral images at room temperature (26 ± 1 ◦C). The system consisted of a computer (Lenovo Tianyi 510 Pro, Lenovo Group Ltd., Beijing, China) installed with data acquisition software (Spectral Image software, Isuzu Optics Corp., Taiwan, China), a spectrograph (ImSpectorV10E, Spectral Imaging Ltd., Oulu, Finland), a 12-bit charged couple device (CCD) camera (HScamera-VIS, Isuzu Optics Corp., Xinzhu, China) with a C-mount lens, an illumination unit of two 150-W tungsten-halogen lamps, and a translation stage (Specim Spectral Imaging Ltd., Oulu, Finland) powered by stepping motor (SC30021A, Zolix Instrument Co, Beijing, China). The spectral solution of the system was 2.8 nm, and the size of incident slit was 30 μm (width) × 14.2 mm (length).

After hyperspectral image acquisition, calibration was conducted using two reference images by the equation as follows:

$$\mathbf{R\_c = (R\_o - D)/(W - D) \times 100\%} \tag{1}$$

where Rc represents calibrated relative reflectance hyperspectral image, Ro denotes acquired original hyperspectral image, D expresses dark reference hyperspectral image with about 0% reflectance, and W stands for white reference hyperspectral image with about 99.9% reflectance. This procedure was directly carried out through the images calibration function within data acquisition software.

#### *2.3. Region of Interests (ROI) Identification and Spectral Extraction*

To isolate pure meat portion in the acquired images, representative regions of interests (ROIs) were first taken based on the calibrated images. ROI was individually determined for each hyperspectral image by applying a corresponding binary mask which was first established using band math procedure. Initially, the low reflectance image at 450 nm was subtracted from the high reflectance image at 890 nm. ROI was formed within the resulting image by thresholding at a constant value of 0.2. The ROIs were completely isolated from backgrounds and edges of the Petri dishes. After that, the mask was built and applied to corresponding hyperspectral images. Then, average spectra were extracted to represent each corresponding sample. All the involved steps were conducted using functions of mask building and ROI generation in ENVI 5.3 (Research Systems Inc., Solutions, Boulder, CO, USA).

#### *2.4. Spectral Pretreatments*

In order to compare and obtain robust and reliable performance, spectral data pretreatment is necessary prior to the development of quantitative models. In this research, a series of pretreatments including normalization, standard normal variate (SNV), multiplicative signal correction (MSC), detrending (detrend), first-order derivative, and second-order derivative (1st and 2nd derivative) were applied in addition to non-preprocessed spectra. Multiplicative interferences of scatter in the spectra can be effectively removed using SNV approach. MSC is a method like SNV which is effectively used in multiplicative variations elimination. Derivatives were usually used to remove baseline

offsets and separate overlapping absorption bands. In this research, derivatives were calculated using second-order polynomial with Savitzky–Golay smoothing by a moving window size of 15 points. Detrending was implemented combined with SNV to suppress the baseline shifting and curvilinearity. Normalization was utilized to present the spectral differences caused by slight optical path variations. The pretreatment or the combinations were all implemented in the Unscrambler X 10.1 (Camo Software Inc., Trondheim, Norway).

#### *2.5. Modeling Method*

PLSR is a reliable linear regression modeling method that has been widely employed in spectral analysis for quantitatively predicting agro-products' quality traits. This method is especially suitable in performing the situation where there is a linear relationship between attributes of the targets and variables. In this study, PLSR models were developed with the dataset in calibration set using a leave-one-out cross-validation (LOOCV) method to prevent data over-fitting. PLSR first projects the spectra onto a few orthogonal factors named latent variables (LVs) [31]. The optimal LVs were determined where lowest root mean square error value of cross-validation was achieved. The PLSR modeling procedure was carried out using the software MATLAB 2013b (MathWorks Inc., Natick, MA, USA) with the PLS toolbox.

#### *2.6. Wavelengths Selection Methods*

Principal component analysis (PCA) is an unsupervised exploratory technique that has been reported to be a powerful tool in dimensionality reduction and multivariate data visualization. The variances of the whole dataset are first explained by PCA, and only a few new orthogonal latent variables that maximize the data variance (called the principal components: PCs) are retained [32]. PCA helps to look for relationships among the samples, and samples in the same class will gather together in PC score plots. If a clear clustering of grouped samples is shown in the PC space by combining two or more PCs, corresponding PC loadings are effective to determine informative wavelengths. Then, pronounced peaks and valleys of the PC loadings are considered to contribute more to the spectral variations of samples with different adulteration levels. The PCA procedure was conducted using the MATLAB 2013b.

Two-dimensional correlation spectroscopy (2D-COS) is a commonly used analytical mathematical formalism. Recently, 2D-COS analysis is highly concerned with identifying a set of spectroscopic data including Raman, visible-infrared, and fluorescence under external perturbation [33]. In terms of generalized 2D-COS, the perturbation can be pressure, concentration, temperature, etc. [34]. To discuss the generated spectrum, synchronous spectrum is diagonal symmetry and there were several autopeaks located at diagonal line. The synchronous spectrum could be used to characterize differences of the spectral intensities at different wavelengths. If the spectral intensities changed sharply with levels at a certain wavelength, there will be a strong autopeak. Thus, in our study, the autopeaks were introduced and utilized as effective wavelengths for identifying the adulteration.

Regression coefficients (RC), which are always used in combination with PLSR modeling method, are also an effective wavelength selection approach. In PLSR models, peaks and valleys at certain wavelengths with dominated RC values indicate a high influence on the response (predicted results) [35]. These spectral variables would be more useful in PLSR modeling and should be chosen for further PLSR models simplification. In our study, wavelengths with high absolute RC values (above the cutoff threshold) from the optimal PLSR model were considered to contribute most in predicting adulteration levels and they were finally adopted.

#### *2.7. Models Performance Assessment*

In order to assess the performance of established models, the following criteria including coefficient of determination in calibration set (Rc 2), cross-validation (Rcv2), and prediction sets (Rp 2), as well as root mean squared error in calibration (RMSEC), cross-validation (RMSECV), and prediction sets (RMSEP) were determined, respectively. Furthermore, residual predictive deviation (RPD) was also evaluated to assess the practical utility of prediction models. If the values of R2 <sup>≥</sup> 0.70 and RPD <sup>≥</sup> 2.00, models were considered to be effective in detecting meat quality and safety [36]. A satisfactory model should perform with results of high values in Rc 2, Rcv2, Rp 2, and RPD as well as low values in RMSEC, RMSECV, and RMSEP.

#### *2.8. Distribution Maps of the Adulterant*

Recognizing adulterate distribution in minced pork is helpful to rapidly observe general adulteration level or if the sample was adulterated. The generation of distribution map is a means of visualization. It is a special advantage of HSI that conventional imaging or NIRS could not achieve. The optimal simplified model can be applied back to predict values in each pixel in the multispectral images at selected wavelengths. After that, the distribution map which pieced all predicted values of pixels together is generated. Therefore, in this research, after the optimal simplified model was confirmed, a colorful image with a linear color scale used for visualizing adulteration levels was displayed. All these steps were performed using a homemade program developed in Matlab 2013b.

#### **3. Results and Discussion**

#### *3.1. Spectral Profiles*

The average reflectance spectra of all the adulterated samples (Figure 1a) and pure minced pork and jowl meat extracted from the corresponding hyperspectral images (Figure 1b) are presented in Figure 1. Different spectra showed similar patterns with certain differences in reflective intensity. As shown, minced pork samples had slightly higher reflective intensity than minced jowl meat samples in spectral region of 400–1000 nm. Although there were few overlaps between spectra of pork and jowl meat, it was still possible to observe certain spectral differences, especially at prominent peaks or valleys. The variations in spectral reflectance among pork and jowl meat were related to the differences in chemical composition, and this implies the adulteration would induce significant alterations to the pure pork samples in a way that can be detected using spectral information.

In the VNIR region, several downwards peaks (absorbance bands) could be observed for pork and jowl meat. The band centered at around 411 nm was associated with the Soret absorption band, which was due to a respiratory pigment of haemoglobin [37]. The valleys of 543 nm and 570 nm can be ascribed to the presentence of deoxymyoglobin and oxymyoglobin, respectively [38], which were responsible for color traits of meat. The 975 nm and 759 nm were attributed to the second and third overtone of O–H stretching mode of water, respectively [39,40]. As for a weak reflectance valley, the independent absorbance band of 842 nm could correspond to the C-H stretching mode of aliphatic compounds [41,42]. The comparison showed that the reflective differences at these bands indicated that minced jowl meat samples had more contents of myoglobin, fat, and water than pure pork. Thus, the identification of jowl meat adulteration in pork is preliminarily concluded to be feasible on the basis of its spectral characteristics difference.

**Figure 1.** Spectral characteristics of prepared samples in the visible and near-infrared region. (**a**) Adulterated samples, (**b**) Mean raw spectra of pure pork and jowl meat.

#### *3.2. PLSR Models*

To determine and quantify the adulteration levels, PLSR models were developed based on raw or pretreated (normalization, SNV, MSC, SNV + Detrend, 1st derivative and 2nd derivative) spectra. A summary of the predictive results is listed in Table 1. As can be seen, spectral data with or without various pretreatments all showed good capability in predicting the adulteration levels. The overall R<sup>2</sup> values were higher than 0.94, RMSE values were lower than 10.2%, and RPD values were higher than 3.1. The small differences observed among RMSEC, RMSECV, and RMSEP values indicated that models were robust and reliable. Overall, HSI coupled with PLSR modeling method provided an innovative way to perform the instant and noncontact prediction for jowl meat adulterated in pork. With regard to a comparison of different pretreatments, the best PLSR modeling results were obtained based on SNV pretreated spectral data with performance of Rp <sup>2</sup> = 0.9549, RMSEP = 7.04%, and RPD = 4.54. Thus, the pretreatment of SNV was adopted to preprocess spectra in subsequent analysis of model simplification and visualization.


**Table 1.** Performance of partial least square regression (PLSR) models for predicting minced jowl meat adulterated in pork.

Notes: PLSR: partial least squares regression; SNV: standard normal variate; MSC: multiplicative scatter correction; LVs: latent variables; Rc 2: coefficient of determination in calibration set; Rcv2: coefficient of determination in cross-validation set; Rp2: coefficient of determination in prediction set; RMSEC: root mean squared error for calibration set; RMSECV: root mean squared error for cross-validation set; RMSEP: root mean squared error for prediction set; RPD: residual predictive deviation.

#### *3.3. Wavelengths Selection*

#### 3.3.1. PCA Explanatory Analysis

Principal component analysis (PCA) is an efficient chemometric method, which provides the interpretation of variances among different data points in spectral analysis. In order to compare and highlight the spectral similarities and differences, PCA was first applied to the whole dataset of 176 spectra. The first three PCs which individually accounted for 80.37%, 9.50%, and 6.26% of the total variance were retained. The reason was that above 95% of the variation could be explained by the first three PCs. Moreover, through trial and error with different PC combinations, PC1 and PC3 were found to be useful in grouping samples into different adulteration levels. Then, the calculated PC scores of data points for different adulteration levels were utilized to create a 2D score plot (Figure 2). In general, data points in the same adulteration level tended to gather together and will be separated from others. The score plot of the combination of PC1 vs. PC3 is shown in Figure 2a. As the adulteration level continued to rise, corresponding samples tended to move along the positive directions of PC1 axis and PC3 axis. However, in low-level adulteration (less than 30%), data clusters were observed to be quite close and overlapped in a certain level. Tracing the root of the above observation, the main chemical composition of homologous pork and jowl meat was too similar so that no distinct separation was displayed if adulteration level was low. In addition, PC2 seemed to mainly express the mutual information of samples with different adulteration levels, which was eliminated in the discrimination.

The PC loading lines of the two effective PCs were analyzed in detail and plotted in Figure 2b. Wavelengths at pronounced peaks and valleys were considered to carry important information in identifying the adulteration and should be selected. As a result, a total of nine wavelengths (440 nm, 491 nm, 545 nm, 560 nm, 570 nm, 632 nm, 686 nm, 752 nm, and 871 nm) were chosen. The valley at 440 nm could be attributed to deoxymyoglobin, the peak at 491 nm is associated with metmyoglobin, and 632 nm could be assigned to sulfmyoglobin [38]. The 676 nm is related to the presentation of redness, and the 871 nm band is relevant with the C-H vibration of hydrocarbons. The wavelengths selected by spectral PCA further confirmed the above results in spectral characteristics that pork and jowl meat were different in color presentation as well as water and hydrocarbon contents.

**Figure 2.** Analysis of effective PC scores and loadings. (**a**) PCA score plot of PC1 vs. PC3, (**b**) Wavelength selection on PC1 and PC3 loading lines.

#### 3.3.2. Two-Dimensional Correction Spectroscopy

The 2D-COS analysis of the obtained 11 average spectra with adulteration levels from 0% to 100% is shown in Figure 3. There were two dominant autopeaks of 491 nm and 632 nm observed at the diagonal line in synchronous spectrum (Figure 3a). Another weak autopeak of 871 nm also occurred which could be clearly seen in the corresponding 3D stereo plot in Figure 3b. The presence of these suggested that intensity at these bands varied seriously with the adulteration levels. Therefore, these three wavelengths were effective in identifying the adulteration levels. It is worth mentioning that these three selected wavelengths were also included in the wavelengths selected by PC loadings. In terms of spectral variables, these three wavelengths are the most important in identifying the adulteration.

**Figure 3.** The 2D-COS spectrum of samples with various adulteration levels. (**a**) Synchronous contour map plot, (**b**) Corresponding synchronous 3D stereo plot.

#### 3.3.3. Regression Coefficients

The regression coefficients (RC) curve from the preferred PLSR model based on SNV pretreated spectra is shown in Figure 4. The cut-off threshold was set to be 5, and only wavelengths at peaks and valleys with higher absolute coefficients than 5 were retained. Finally, a total of 10 (433 nm, 450 nm, 481 nm, 558 nm, 578 nm, 594 nm, 634 nm, 661 nm, 889 nm, and 948 nm) discontinuous wavelengths were deemed as the most effective wavelengths in PLSR models' development for quantitatively predicting jowl meat adulteration in minced pork. These wavelengths were different from the above ones but also reasonable due to the consideration of targeted prediction values.

Based on the above three wavelengths selection methods, the number of variables reduced significantly by at least 96.5% to at most 98.9%. The retained wavelengths could be utilized in developing a robust model and further a low-cost multispectral imaging system. Therefore, all the three groups of effective wavelengths could be the basis for comparison in developing the simplified PLSR models.

**Figure 4.** Regression coefficients of the optimal PLSR model based on full spectra.

#### *3.4. Multispectral Models Development*

In order to further eliminate the useless wavelengths and optimize the processing time in computing, simplified PLSR models based on selected wavelengths were established. Wavelengths selected by 2D-COS, PC loadings, and RC methods were individually set as inputs of the simplified PLSR models, and the overall results are displayed in Table 2. As can be seen, the results obtained using selected wavelengths slightly decreased compared with the full spectra. This phenomenon indicated that selected wavelengths were effective and the eliminated variables also contained little information in determining the adulteration. Compared with 2D-COS-PLSR and PC loadings-PLSR models, RC-PLSR model achieved the best performance with Rp <sup>2</sup> of 0.9063, RMSEP of 13.93%, and RPD of 2.30. It indicated that the 10 wavelengths selected by RC method were the most critical in identifying jowl meat adulteration in pork. On the contrary, the three wavelengths selected by 2D-COS performed not that well mainly because that they were not informative enough. As an extension, the additional six more wavelengths selected using PC loadings significantly improved the prediction accuracy by showing the Rp <sup>2</sup> of 0.7475, RMSEP of 17.31%, and RPD of 1.85. What is more, 2D-COS and PC loadings selected wavelengths using X-variables (spectra) only, while RC was based on PLSR model which decomposed both X-and Y-variables (adulteration levels) in the LVs calculation. RC built the optimal relationship between spectral data and adulteration levels compared with 2D-COS and PC loadings. Therefore, the multispectral RC-PLSR model was finally chosen for further visualization steps.


**Table 2.** Performance of simplified PLSR models based on wavelengths selected by three methods.

*Notes:* PLSR: partial least squares regression; 2D-COS: two-dimensional correction spectroscopy; PC: principal component; RC: regression coefficients; LVs: latent variables; Rc 2: coefficient of determination in calibration set; Rcv2: coefficient of determination in cross-validation set; Rp2: coefficient of determination in prediction set; RMSEC: root mean squared error for calibration set; RMSECV: root mean squared error for cross-validation set; RMSEP: root mean squared error for prediction set; RPD: residual predictive deviation.

The determination of the limit of detection (LOD) is an important step which investigated if the lowest adulteration concentration can be detected with the HSI methodology. The LOD was commonly calculated to evaluate the sensitivity of detection methods by the following equation [43].

$$\text{LOD} = 2\delta\_\text{b/S} \tag{2}$$

where δ<sup>b</sup> indicates the standard deviation (SD) of the background response and *S* denotes the sensitivity by the ratio of the predicted adulteration levels to the reference values (namely the slope of the calibration line).

The performance of the preferred multispectral RC-PLSR model with error bar is illustrated in Figure 5. The results of this optimal model were visualized, and it could be seen that all the adulterated samples could be detected (predicted values were above 0%) whether in the calibration or prediction set. The LOD in independent prediction set calculated by the above Equation (2) was found to be 6.50%. However, the aim of meat adulteration is generally to make profit so that adulteration is always performed to be more than 10% [44]. Thus, the LOD in this research proved that it was competent in detecting jowl meat adulteration in pork by the HSI system.

**Figure 5.** The performance of preferred multispectral RC-PLSR model with error bar in (**a**) Calibration set, and (**b**) Prediction set.

#### *3.5. Visualization of the Adulteration Levels*

As known, each sample was represented by the average spectra extracted from corresponding ROI, and the targeted adulteration level was only indicated by one value. However, there was abundant spatially distributed information in hyperspectral images. In this research, the adulteration level at each pixel in one hyperspectral image was predicted so that distribution was visualized for a quick view. The optimal RC-PLSR model was first applied to the multispectral images recombined at selected

wavelengths. The false color images (left column) and corresponding predicted distribution maps (right column) of the samples with different adulteration levels in prediction set are shown in Figure 6. The false color image in Figure 6a was composited by setting the images at 700.9 nm, 545.4 nm, and 436.4 nm as R (red), G (green), and B (blue) channels using ENVI software. They were quite close to the true color image at the primary RGB colors' wavelengths (700 nm, 546.1 nm and 435.8 nm). As can be seen, actual adulteration level is difficult to be recognized by naked eyes in Figure 6a. Distribution maps expressed how the adulteration varied from sample to sample and even from pixel to pixel within one sample and were generated to be shown in Figure 6b. The linear color bar located in the right side from black to red indicated different adulteration levels from 0% to 100% accordingly. There was a clear tendency of color gradient with the increasing adulteration level so that adulteration was easily distinguishable in distribution maps.

**Figure 6.** Distribution maps of jowl meat adulteration in pork at pixel level. (**a**) False color images, (**b**) Distribution maps.

In further steps, the reliability of visualization was verified through generating distribution maps for the two control samples. Two false color images were in the up row and corresponding distribution maps were displayed in the down row in Figure 7. For control sample one, the results of the distribution map showed that four fan-shaped areas were in different colors. The upper left area was mainly in red, and the lower left showed red and yellow. The upper right exhibited a wide range of colors while the lower right presented general blue. With regard to control sample two, the distribution map gave a display of half black and half yellow. All the prediction maps were generally consistent with the actual situation. The effectiveness of the established RC-PLSR model and visualization procedure was thus proved again.

**Figure 7.** Visualization for control samples.

#### **4. Conclusions**

This study was motivated by the requirement for rapidly and nondestructively identifying one common adulterant of minced jowl meat in minced pork. Our attempts explored the application of HSI combined with chemometrics and wavelength selection algorithm to quantify and visualize the adulteration. In particular, simplified RC-PLSR model developed by 10 key wavelengths gave the best performance; LOD achieved 6.50%. The visualization of adulteration levels was successfully performed based on RC-PLSR model. Predicted colorful distribution maps were generated to make it convenient in observation. To test the validity of visualization, two known distributed samples were predicted and expected corresponding maps were displayed. The overall results suggested that HSI had the potential to identify minced jowl meat adulteration in pork without any prior physical or chemical analysis information. This technique could provide more detailed visualization information than conventional imaging and NIRS in identifying adulteration levels. However, more studies related to a large number of samples and different breeds in sampling should be conducted. In addition, further work will focus on identifying a variety of commonly used adulterates in pork based on a few wavelengths to achieve the goal for portable application.

**Author Contributions:** H.J. conceived and designed the experiments. H.J., F.C., and M.S. conducted experiments and analyzed the data. H.J. wrote the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Scientific Research Foundation for Advanced Talents of Nanjing Forestry University.

**Acknowledgments:** The authors would like to thank organizations or individuals who had contributed to this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
