**Author Contributions:** Conceptualization, methodology, writing—original draft preparation, A.H.; writing original draft preparation, I.M.; writing—original draft preparation, revision, W.F.S.; writing—original draft preparation, H.T.T.; writing—original draft preparation, L.L.; writing—original draft preparation, H.-Y.K.; project administration, supervision, manuscript revision, H.N.; writing—original draft preparation, A.B.; writing—original draft preparation, A.A.-K.; writing—original draft preparation, M.S.; writing—original draft preparation, E.S.; writing—original draft preparation, S.G.; writing—original draft preparation, D.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This work was supported by the Norwegian Institute of Food, Fisheries, and Aquaculture Research (Nofima) through a Strategic Research Initiative (Spectec Project): Rapid and Non-Destructive Measurements to Enable Process Optimization.

**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*

## **Flash Gas Chromatography in Tandem with Chemometrics: A Rapid Screening Tool for Quality**

#### **Sara Barbieri 1 , Chiara Cevoli 1 , Alessandra Bendini 1,\* , Beatriz Quintanilla-Casas 2,3 , Diego Luis García-González 4 and Tullia Gallina Toschi 1**

**Grades of Virgin Olive Oils**

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +39-0547-338121

Received: 11 June 2020; Accepted: 29 June 2020; Published: 2 July 2020

**Abstract:** This research aims to develop a classification model based on untargeted elaboration of volatile fraction fingerprints of virgin olive oils (*n* = 331) analyzed by flash gas chromatography to predict the commercial category of samples (extra virgin olive oil, EVOO; virgin olive oil, VOO; lampante olive oil, LOO). The raw data related to volatile profiles were considered as independent variables, while the quality grades provided by sensory assessment were defined as a reference parameter. This data matrix was elaborated using the linear technique partial least squares-discriminant analysis (PLS-DA), applying, in sequence, two sequential classification models with two categories (EVOO vs. no-EVOO followed by VOO vs. LOO and LOO vs. no-LOO followed by VOO vs. EVOO). The results from this large set of samples provide satisfactory percentages of correctly classified samples, ranging from 72% to 85%, in external validation. This confirms the reliability of this approach in rapid screening of quality grades and that it represents a valid solution for supporting sensory panels, increasing the efficiency of the controls, and also applicable to the industrial sector.

**Keywords:** virgin olive oil; quality; volatile compounds; sensory analysis; chemometrics

### **1. Introduction**

The official methodology for sensory evaluation of virgin olive oils (VOOs), known as a panel test, is a fundamental tool to assess the quality of products that cannot be replaced by instrumental methods, considering that the overall and complex perceptual attributes (e.g., fruity and defects) are the indicators of the quality of VOOs. Despite its proven effectiveness in evaluating the quality grades of samples, tested in EU countries since 1991 [1,2], the scientific community has highlighted some drawbacks on its application that are mainly related to the following: (i) the reproducibility of results among different panels; (ii) critical attribution of the category when, e.g., a defect is borderline; (iii) costs, assessor fatigue and other limitations associated with a method working with humans.

Specifically, according to decisions taken at International Olive Council (IOC) level, the Reg. (EU) 1348/2013 [3] recommends the number of oils to be assessed by the sensory panels, fixing a maximum number of four samples at each session. Moreover, a maximum of three sessions per day is specified, to leave enough time between a session and another, thus avoiding the contrast effect that could be produced by immediately tasting sequences of samples. These specifications strongly limit the number of samples that can be assessed by one panel per day. On the other hand, to enhance panel skills in recognizing, identifying, and quantifying sensory attributes, the introduction of new artificial reference materials (obtained by chemical or biotechnological approaches), could improve the proficiency of the individual panels and their global alignment by overcoming some limitations associated with a natural matrix (e.g., limited amounts available, difficultly obtaining, low homogeneity year by year) and offering advantages such as preparation in each laboratory, reproducibility over time, possibility of purchase, and therefore their availability for the market.

In this context, the development of an instrumental method for rapid screening of quality grades of samples (extra virgin olive oil, EVOO; virgin olive oil, VOO; lampante olive oil, LOO) could represent a solution to support sensory panels (particularly for large private industries), decreasing their daily work by reducing the samples that need to be assessed (e.g., by excluding those definitely compliant), with a consequent increase in the efficiency of quality controls and reducing the number of samples that need to be controlled.

In this way, improvement of the activity of sensory panels, whose work remains central to ensuring the quality of the product, would be achieved by focusing sensory analysis only on uncertain samples (i.e., borderline oils between two product categories that can be the object of disagreement among panels).

It is well known that volatile compounds are crucial to determine VOO quality and that they are responsible for the different VOO sensory profiles [4–6]; their determination in a rapid way (e.g., screening method) could support sensory analysis and represents one of the current challenges in the olive oil sector where fast, accurate, and easy-to-use approaches providing real-time results are required.

Recently, different analytical techniques combined with chemometric statistical approaches have been proposed to predict sensory information [7–9].

Alongside the traditional techniques (targeted) in which specific and selected molecular markers are monitored during the analysis to assess the presence or absence of compounds and their quantification, untargeted analyses, based on a holistic approach and able to provide information such as a spectral fingerprint, giving a simplified and overall picture of the food under analysis, have gained an increasing relevance over the last years [10].

Among the latter, different analytical methods for determination of volatile compounds combined with multivariate chemometric techniques for VOOs quality testing have been described in the literature and proposed to the industrial sector as fast and high throughput screening techniques [9,11–18].

In particular, as an alternative to headspace gas chromatography-mass spectrometry (HS-GC/MS), which is the most widely used technique to quantify and characterize the profiles of volatile compounds of VOOs thanks to its high sensitivity and selectivity, the application of the HS-GC ion mobility spectrometry (HS-GC-IMS) has been proposed. This technique combines high selectivity and sensitivity with high robustness and cost-efficiency, and has given promising results in discriminating VOOs according to quality grades [9,11,12,14,18] or geographical origin [13,15].

The need to support organoleptic analysis was also reported in a specific call of the Horizon 2020 EU program (H2020-SFS-14a-2014) and is one of the main objectives of the OLEUM project (Horizon 2020, Grant Agreement No. 635690). In the framework of this project, two analytical instrumental techniques, headspace-solid phase micro extraction–gas-chromatography/mass spectrometry (HS-SPME-GC/MS) and flash gas chromatography (FGC) based on the determination of volatile compounds, have been proposed as the most promising rapid screening methods that can support sensory panels in the determination of quality grades.

In a recent work by Quintanilla-Casas and co-authors (2020) [17], the results obtained with HS-SPME-GC/MS with a fingerprinting approach to classify VOO categories has been demonstrated. Herein, a classification model based on minor fraction fingerprints that is able to predict the commercial

category of olive oil samples (EVOOs, VOOs, LOOs) obtained by FGC is presented. The FGC is an innovative analytical approach for analysis of volatile compounds of VOOs based on the FGC separation: the headspace of VOOs, previously conditioned, is sampled by a syringe, the volatile organic compounds are adsorbed on a Tenax trap and subsequently desorbed by rapid heating, and, finally, transferred to a FGC step. The elution of analytes runs in parallel using two metal capillary columns with different polarity of the stationary phase. This gives rise to slight differences in the separation capability of molecules that are detected by a flame ionization detector (FID) located at the end of each column.

The main advantage of the FGC technique is its short analysis time (total separation time is 100 s); moreover, its application associated with sensory analysis for calibration and chemometric tools is promising to support the work of panel tests in discriminating samples of different product categories. A classification model, once built, could be easily applied in any laboratory or industry.

The effectiveness of this technique is already demonstrated by previous works aimed to differentiate VOOs according to their geographical origin declared by labels such as "100% Italian" and "non-100% Italian" oils [19] or "EU" and "extra-EU" [16].

The aim of this study was to classify VOOs according to quality grade, combining FGC data with the multivariate classification technique partial least squares discriminant analysis (PLS-DA). To provide robustness to our model, a set of 331 oils belonging to the three different commercial categories (EVOO, VOO, LOO) involving two harvesting/production years was analyzed. The adopted validation protocol (repeatability and reproducibility tests) and related performance are also shown.

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

#### *2.1. Olive Oil Samples*

An initial set of 334 EVOOs, VOOs, and LOOs oils representative of the most common olive cultivars, geographical origin, sensory positive attributes, and sensory defects were sampled. Specifically, in addition to a first set of 180 oils collected during the first year of the OLEUM project (2016–2017 olive season), another set of 154 samples (2017–2018 olive season) was collected and analyzed during the second year (Tables S1–S4 in the Supplementary Materials).

The panel test method was carried out by six panels involved in the OLEUM project as described by Barbieri et al. 2020 [20] and sensory data were expressed as mean of medians. The procedure deals with possible disagreement between panels with a decision tree in order to have definitive classification of samples in which definitive agreement is reached. In agreement with the sensory results reported in Tables S1 and S2 (Supplementary Materials), in the first year of the project 178 of 180 samples were immediately classified by panels (54 EVOO, 78 VOO, and 48 LOO). Classification was not possible for only two samples (UN\_10, UP\_14), as agreement among panels was not reached on the category (V/L). The sensory evaluation of oils from the second sampling allowed classification of 153 oils (69 EVOO, 51 VOO and 33 LOO); 1 sample was not classified due to an anomalous lemon smell (ZRS\_1) and was therefore excluded from the set [20]. For these reasons, the classification model was built on 331 samples.

The oils collected were representative of possible commercial samples and borderline samples that can be the object of disagreement between panels in terms of sensory characteristics. Different aliquots of the samples, stored in the lab at 10–12 ◦C (for sensory analysis) and at −18 ◦C (for instrumental analysis), were reconditioned at room temperature before analysis.

#### *2.2. Analytical Conditions*

The FGC system (FGC-E-nose Heracles II, AlphaMos, Toulouse, France) is based on the technology of ultra-fast gas-chromatography.

The FGC is equipped with two columns working in parallel: a non-polar column (MXT5: 5% diphenyl, 95% methylpolysiloxane, 10 m length and 180 µm diameter) and a polar column (MXT-1701: 14% cyanopropylphenyl/86% dimethyl polysiloxane, 10 m length, 180 µm diameter). At the end of each column, a FID detector is placed and the acquired signal is digitalized every 0.01 s.

The analytical conditions applied were the same described by Melucci et al. 2016 [19]. The only difference was related to the temperature of the conditioning step of the samples before injection: the vial is placed in the auto-sampler (HS 100, CTC Analytics), which moves it in a shaker oven where it remains for 20 min at 40 ◦C, shaken at 500 rpm.

#### *2.3. Validation Protocol*

To confirm that the analytical procedure employed has performance capabilities consistent with the required application, a validation strategy for non-targeted approaches was performed.

A QC (quality control) sample, representative of the qualitative and quantitative VOO volatile composition (presence of volatile compounds along the entire interval of the chromatogram), was used. In this study, the QC sample was obtained by pooling the same volume of three case-control samples (1 EVOO, 1 VOO with median of 1.9 for fusty-muddy defect, and 1 VOO with a median of 2.5 for rancid defect) and seven replicates were taken into consideration.

The quality of the instrumental performance intended for fingerprinting analysis was checked by the calculation of the relative standard deviation (RSD) as proposed by the Food and Drug Administration [21]. Specifically, the repeatability (intra-day repeatability and inter-day repeatability performed according to EC 657/2002) [22] of the chromatographic signal evaluated in terms of RSD% of each chromatogram data point, with intensities above noise signal of the replicates of the same QC samples, was considered [23,24].

Prior to RSD calculation, data were aligned using the COW algorithm (correlation optimized warping) [25] and autoscaled (mean-centering followed by division of variable by the standard deviation of that column) to correct shifts in retention time and possible differences in the signal amplification of the instrument. All elaborations were made using PLS Toolbox for Matlab (MatlabR2018a®) (Natick, MA, USA).). For calculation of RSD% for each chromatogram data point, the evaluation and exclusion of noise signal is carried out to avoid considering non-relevant RSD%.

For precision, the FDA recommends a RSD not higher than 15% regarding the analytical variability for target analysis, except for concentrations close to the detection limit where a RSD of 20% is acceptable (FDA Bioanalytical Method Validation-Guidance for Industry, 2018). This, in agreement with the trend described by the Horwitz equation for targeted methods [26], demonstrates that the repeatability is strongly correlated with the intensity of the variables.

Although fingerprinting represents a different analytical approach and more variation is expected when doing untargeted analysis, these guidelines are used as a benchmark towards repeatability evaluation. Specifically, for intra-day repeatability, the acceptance criteria were as follows: more than 90% of signals with RSD < 15%; more than 95% of signals with RSD < 20% and distribution of RSD% vs. signal intensity in accordance with the Horwitz equation. For inter-day repeatability or within-lab reproducibility, the acceptance criteria were as follows: more than 85% of signals with RSD < 15%, more than 90% of signals with RSD < 20% and distribution of RSD% vs. signal intensity in accordance with Horwitz's equation.

In addition, the examination of system performance by checking the signal to noise ratio in standard solutions (instead of the evaluation of representative VOO profiles) to facilitate the assessment and comparison of method sensitivity for other laboratories was proposed. The sensitivity of the analytical system was evaluated by analyzing 2 g of each standard solution in refined olive oil (ethanol 0.05 mg·kg−1 , CAS Number 64-17-5; assay ≤ 97.2%; density 0.789 g/mL at 25 ◦C; hexanal, 0.1 mg·kg−1 CAS Number 66-25-1; assay ≥ 95% (GC); density 0.815 g/mL at 25 ◦C; (*E*)-2-hexenal, 0.75 mg·kg−1 CAS Number 6728-26-3; assay ≥ 97.0% (GC); density 0.846 g/mL at 25 ◦C). The S/N (S = intensity of the peak of the compound; N = mean intensity of the noise measured considering the baseline of the chromatographic zone between 43 and 50 s) for the selected analytes in the chromatograms should be >3 (acceptance criteria).

#### *2.4. Classification Models*

In order to predict the assignment of samples to a specific quality grade, full chromatograms were used to develop classification models. The raw data of each chromatogram, for a total of 19,900 points, were aligned by the COW algorithm and autoscaled using PLS Toolbox for Matlab (MatlabR2018a®). Subsequently, the noise was excluded and 8401 points were consecutively selected from first to last peak observed in the chromatogram.

Subsequently, PLS-DA (partial least square discriminant analysis) models [27] were built by using the intensity values of the points as variables X (matrix X), while the commercial categories (EVOO, VOO, LOO) were considered as variable Y. In particular, classification models with 2 categories were developed in sequence: EVOO vs. no-EVOO followed by VOO vs. LOO and LOO vs. no-LOO followed by EVOO vs. VOO, as proposed by Quintanilla-Casas et al. 2020 [17].

The sample dataset was split in calibration (venetian blinds cross validation, including 75% of the samples) and external validation set (25% of the samples) by using the Kennard–Stone method [28]. The dataset was deposited for possible consultation in an on-line repository [29].

The threshold value able to identify the belonging category of each sample into one of the groups was defined by using a probabilistic approach based on Bayes's rule [30]. Finally, to assess the goodness of the method, the receiver operating characteristic (ROC) curves were evaluated.

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

#### *3.1. Performance of FGC*

Most of the procedures proposed in the literature for validation of non-targeted methods focus on post-analytical data treatment and validation of statistical models. Nevertheless, a few studies have investigated control procedures as well as performance criteria and requirements to ensure the consistence of the analytical signal (fingerprint) [24,31].

Conventional performance criteria adopted for targeted methods are not applicable as such to fingerprinting methods. Fingerprinting methods intended for sample classification are not aimed at identification and quantification of analytes, but on finding distinctive patterns that are specific for a given food category (i.e., VOO commercial category) in raw analytical signals (i.e., chromatograms). Therefore, the main constraint of the fingerprinting analytical method is to provide a repeatable and reproducible signal with sufficient sensitivity to collect the information from samples for the final purpose of the method, i.e., quality classification.

For evaluation of intra-day repeatability, the pooled QC sample was analyzed by the same operator with the same equipment and in the same instrument operative conditions within the same day. For each variable (data points), mean value, SD, and RSD% were calculated considering the seven replicates. More than 97.5% of signals presented RSD < 10%, while it achieves 99.8% in correspondence of RSD < 20% (Table 1). To analyze the variability as related to the magnitude of the variables, RSD% was plotted versus signal intensity (data not shown). As expected, data points with RSD > 10% are characterized by low values of intensity. This is in agreement with the trend described by the Horwitz equation for targeted methods [26].

In the case of the inter-day repeatability (within-lab reproducibility), seven replicates of the pooled QC sample were analyzed by the same operator with the same equipment but on different days, consequently involving different environmental conditions, and the mean value, SD, and RSD% were calculated. More than 91% and 99.4% of the signals presented RSD < 10% and RSD < 20%, respectively (Table 1). A relation between intensity and RDS% was also observed in this study, similarly to that previously observed in the intra-day repeatability test.

As the fingerprinting approach intended for sample classification is not aimed in determining the concentration of single analytes, limits of detection or quantification cannot be calculated for the analytical outcome. However, the analytical method needs to be sufficiently sensitive to allow detection of minor constituents to avoid missing any valuable information.

**Table 1.** Frequency of each relative standard deviation percentage (RSD%) class obtained for intra-day and inter-day repeatability evaluated on the quality control (QC) sample.

On this basis, the method's sensitivity needs to be set as a reference parameter to be evaluated in the validation process. A target-type strategy applied to standard solutions was proposed.

Standard solution compounds were chosen as most representative of the qualitative and quantitative volatile composition of VOOs, especially regarding the presence of volatile compounds over the entire interval of the chromatogram considered in fingerprinting analysis. Differences between the concentrations used for each compound are related to their different amounts generally present in a VOO sample. Results of the S/N are reported in Table 2.

**Table 2.** Concentration (mg·kg−1 ) of each compound included in the standard solution used for method's sensitivity evaluation and related S/N. The standard mix were prepared by spiking refined olive oil with each compound and analysed by flash gas chromatography (FGC). S = intensity of the peak of the compound; N = mean intensity of the noise measured considering the baseline of the chromatographic zone between 43 and 50 s.

#### *3.2. Classification Models*

A fingerprinting approach involving chemometric elaboration of the entire profiles in volatile molecules without identification and quantification was applied.

Two different classification strategies were taken into account: (i) a classification model able to discriminate EVOO and no-EVOO samples, followed by a model to classify VOO vs. LOO samples; (ii) a classification model able to discriminate LOO and no-LOO samples, followed by a model to classify VOO vs. EVOO samples.

The results, in terms of percentage and number of correctly classified samples, are reported in Table 3 for cross and external validation, respectively. Regarding the first classification strategy, the percentages of correctly classified samples ranged from 72 to 89% and from 72 to 85%, for cross and external validation, respectively. In particular, the best results were obtained during the second step useful to discriminate VOO vs. LOO. For the second strategy, conceptually more correct in terms of sequence because it first discriminates LOO which are not edible if not refined, the percentage ranged from 78 to 92% and from 73 to 85%, for cross and external validation, respectively. In this case, the highest percentages were reached using the first PLS-DA model (LOO vs. no-LOO). Furthermore, this latter model was the best of all PLS-DA models developed.

oil, LOO = lampante olive oil.

In general, the percentages are in the same range as those obtained by other authors who proposed chemometric models to discriminate VOO quality grades according to their volatile profile analyzed by different instrumental techniques [9,17]. In general, the percentages are in the same range as those obtained by other authors who proposed chemometric models to discriminate VOO quality grades according to their volatile profile analyzed by different instrumental techniques [9,17].

TOTAL: 129/156 = 83% TOTAL: 41/51 = 78% TOTAL: 158/189 = 84% TOTAL: 47/60 = 78%

*Foods* **2020**, *9*, x FOR PEER REVIEW 8 of 12

**Table 3.** Results in terms of percentage and number of samples correctly classified in cross and external validation of the two classification strategies applied based on the partial least squaresdiscriminant analysis (PLS-DA) sequential model. EVOO = extra virgin olive oil; VOO = virgin olive

**1st CLASSIFICATION STRATEGY 2nd CLASSIFICATION STRATEGY 1st Step: EVOO vs. no-EVOO 1st Step: LOO vs. no-LOO**  Cross validation External validation Cross validation External validation EVOO: 70/90 (78%) EVOO: 26/32 (81%) LOO: 50/61 (81%) LOO: 17/20 (85%) No-EVOO: 132/164 (81%) No-EVOO: 37/48 (77%) No-LOO: 172/188 (92%) No-LOO: 55/65 (85%) TOTAL: 202/254 = 80% TOTAL: 63/80 = 79% TOTAL: 222/249 = 89% TOTAL: 72/85 = 85% **2nd Step: VOO vs. LOO 2nd Step: VOO vs. EVOO**  Cross validation External validation Cross validation External validation VOO: 88/99 (89%) VOO: 22/26 (85%) VOO: 84/95 (88%) VOO: 23/27 (85%)

The ROC curves (Figure 1) evaluated the sensitivity (number of samples predicted as in the class divided by number actually in the class) and the specificity (number of samples predicted as not in the class divided by actual number not in the class) of all PLS-DA models (external validation) [16]. In particular, the area under the curve (AUC) identifies the degree of discrimination (ranged 0.8148 to 0.8899) and suggests that all the models are characterized by a good degree of discrimination. The ROC curves (Figure 1) evaluated the sensitivity (number of samples predicted as in the class divided by number actually in the class) and the specificity (number of samples predicted as not in the class divided by actual number not in the class) of all PLS-DA models (external validation) [16]. In particular, the area under the curve (AUC) identifies the degree of discrimination (ranged 0.8148 to 0.8899) and suggests that all the models are characterized by a good degree of discrimination.

**Figure 1.** Receiver operating characteristic (ROC) curves of all developed PLS-DA models used to discriminate samples according to quality grade; the red circles identify the sensitivity (number of samples predicted as in the class divided by number actually in the class) and the specificity (number of samples predicted as not in the class divided by actual number not in the class) of the models. EVOO = extra virgin olive oil; VOO = virgin olive oil, LOO = lampante Olive Oil. **Figure 1.** Receiver operating characteristic (ROC) curves of all developed PLS-DA models used to discriminate samples according to quality grade; the red circles identify the sensitivity (number of samples predicted as in the class divided by number actually in the class) and the specificity (number of samples predicted as not in the class divided by actual number not in the class) of the models. EVOO = extra virgin olive oil; VOO = virgin olive oil, LOO = lampante Olive Oil.

The results of all the models (cross and external validation), in term of probability of belonging to the correct class, are shown in Figure 2. The threshold value was fixed at 0.5, corresponding to a probability of 50%: a sample classified with a probability lower than this is considered as not correctly grouped [32].

The definition of a probability level, ranging from 50% to 100%, could be a means of identifying uncertain samples that need to be checked by sensory evaluation. In other words, the samples classified with a probability lower than the selected probability level should be submitted to panel test. These procedures would reduce the amount of the samples analyzed by the panel, but at the same time, it would insure the accuracy of the classification.

same time, it would insure the accuracy of the classification.

grouped [32].

The results of all the models (cross and external validation), in term of probability of belonging to the correct class, are shown in Figure 2. The threshold value was fixed at 0.5, corresponding to a probability of 50%: a sample classified with a probability lower than this is considered as not correctly

The definition of a probability level, ranging from 50% to 100%, could be a means of identifying uncertain samples that need to be checked by sensory evaluation. In other words, the samples classified with a probability lower than the selected probability level should be submitted to panel

**Figure 2.** Class prediction probability of all samples used to develop the models, in cross and external validation (grey area). Step 1—EVOO (green star) vs. no-EVOO (blue square); step 2—VOO (yellow diamond) vs. LOO (red circle); step 1—LOO (red circle) vs. no-LOO (yellow square); step 2—VOO (yellow diamond) vs. EVOO (green star). EVOO = extra virgin olive oil; VOO = virgin olive oil, LOO = lampante olive oil. **Figure 2.** Class prediction probability of all samples used to develop the models, in cross and external validation (grey area). Step 1—EVOO (green star) vs. no-EVOO (blue square); step 2—VOO (yellow diamond) vs. LOO (red circle); step 1—LOO (red circle) vs. no-LOO (yellow square); step 2—VOO (yellow diamond) vs. EVOO (green star). EVOO = extra virgin olive oil; VOO = virgin olive oil, LOO = lampante olive oil.

**4. Conclusions**  Despite the undisputed validity of the panel test, its application is time consuming and expensive. Accordingly, companies and private and public quality control labs could benefit from **Table 3.** Results in terms of percentage and number of samples correctly classified in cross and external validation of the two classification strategies applied based on the partial least squares-discriminant analysis (PLS-DA) sequential model. EVOO = extra virgin olive oil; VOO = virgin olive oil, LOO = lampante olive oil.

#### **4. Conclusions**

Despite the undisputed validity of the panel test, its application is time consuming and expensive. Accordingly, companies and private and public quality control labs could benefit from robust instrumental pre-classifications, which would reduce the number of samples that have to be assessed by panels, or at least prioritize their assessment.

For this reason, the development of rapid screening methods to support the official panel test, to analyze olive oils and differentiate their quality grades, is one of the challenges in the olive oil sector, as reported in the EU framework program Horizon 2020.

In this work, FGC combined with the multivariate statistical technique was applied to discriminate samples according to different quality grades (EVOO, VOO and LOO; examples of GC traces for EVOOs and LOOs are shown in Figure S1 of the Supplementary Materials). The analytical technique proposed herein for fingerprinting olive oils combined with chemometrics was effective in reducing data complexity and time to obtain a response; this rapid screening tool could be adopted for a quick pre-classification of the quality grades, e.g., by control laboratories in companies of the OO sector, before buying or blending EVOOs.

In order to propose a robust chemometric model, a large set of samples (*n* = 331) involving two different harvesting/production years, the most common olive cultivars, geographical origin, sensory positive attributes, and sensory defects, was analyzed. In addition, a validation protocol was adopted for evaluate the reliability of the results.

The proposed analytical fingerprinting method provided repeatable and reproducible signals with sufficient sensitivity to collect valuable information about samples.

FGC associated with the two-category sequential classification model is promising to support sensory analysis in discriminating samples of different product categories. Among the proposed classification strategy, the second (1st step: LOO vs. no-LOO; 2nd step: VOO vs. EVOO) was the best of all PLS-DA models developed with percentages of correctly classified samples ranging from 78 to 92% and from 73 to 85%, for cross and external validation, respectively.

This analytical approach is very fast, and, in fact, only around 200 s are needed to analyze a single sample. The classification model, built by using a high number of robust samples classified by sensorial analysis and representative of the commercial variability (here we used a decision tree and six panels to ensure their classification) is easily applicable in any laboratory or industry.

Future studies could be addressed to the implementation of this methodology, even in relation to an increasing interest of the food sector towards volatile compounds and more widespread use of instruments such as FGC, which are less common in quality control laboratories. An even wider sampling phase including other variables among oils, since they are natural products, could lead to a better control of classifications and would lead to implementation of this technique to a broader extent. Lastly, the use of other statistical approaches, such as nonlinear techniques, could be investigated in order to improve the results of classification.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2304-8158/9/7/862/s1, Table S1: Sensory results of samples from the first year. Table S2: Sensory results of samples from the second year. Table S3: available information on samples collected and evaluated during the first year of the Oleum project. Table S4: available information on samples collected and evaluated during the second year of the Oleum project. Figure S1: overlapping of the GC traces of extra virgin (EVOO) and lampante (LOO) samples.

**Author Contributions:** Conceptualization, S.B., A.B. and T.G.T.; Formal analysis, S.B. and C.C.; Data curation, S.B. and C.C.; Writing—original draft preparation, S.B.; Writing—review and editing, S.B., C.C., A.B., B.Q.-C.; Supervision, T.G.T. and D.L.G.-G.; Funding acquisition, T.G.T. All authors have read and agree to the published version of the manuscript.

**Funding:** This work is supported by the Horizon 2020 European Research project OLEUM "Advanced solutions for assuring the authenticity and quality of olive oil at a global scale", which received funding from the European Commission within the Horizon 2020 Programme (2014–2020), grant agreement No. 635690.

**Acknowledgments:** The information expressed in this article reflects the authors' views; the European Commission is not liable for the information contained herein. We are grateful to all producers who provided us with VOOs for this study as well as the panel members who performed sensory analysis of VOOs from each institution involved: Eurofins Analytik GmbH, Hamburg, Germany; Institute of Agriculture and Tourism, Porˇec, Croatia; Institut des Corps Gras, Pessac, France; Alma Mater Studiorum-Università di Bologna; Science and Research Centre Koper, Slovenia and Ulusal Zeytin ve Zeytinyăgı Konseyi, Izmir, Turkey.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

## **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/).

## **Genetic Identification of the Wild Form of Olive (***Olea europaea var. sylvestris***) Using Allele-Specific Real-Time PCR**

### **Christina I. Kyriakopoulou and Despina P. Kalogianni \***

Department of Chemistry, University of Patras, 26504 Rio Patras, Greece<EMAIL_ADDRESS>**\*** Correspondence<EMAIL_ADDRESS>
Received: 28 February 2020; Accepted: 7 April 2020; Published: 9 April 2020

**Abstract:** The wild-type of olive tree, *Olea europaea var Sylvestris* or oleaster, is the ancestor of the cultivated olive tree. Wild-type olive oil is considered to be more nutritious with increased antioxidant activity compared to the common cultivated type (*Olea europaea* L. *var Europaea*). This has led to the wild-type of olive oil having a much higher financial value. Thus, wild olive oil is one of the most susceptible agricultural food products to adulteration with other olive oils of lower nutritional and economical value. As cultivated and wild-type olives have similar phenotypes, there is a need to establish analytical methods to distinguish the two plant species. In this work, a new method has been developed which is able to distinguish *Olea europaea var Sylvestris* (wild-type olive) from *Olea europaea* L. *var Europaea* (cultivated olive). The method is based, for the first time, on the genotyping, by allele-specific, real-time PCR, of a single nucleotide polymorphism (SNP) present in the two olives' chloroplastic genomes. With the proposed method, we were able to detect as little as 1% content of the wild-type olive in binary DNA mixtures of the two olive species.

**Keywords:** *Olea europaea var Sylvestris*; oleaster; olive; olive oil; real-time PCR; adulteration; SNP; DNA

### **1. Introduction**

The wild form of the olive tree, formally named *Olea europaea var Sylvestris* or oleaster, is considered to be one of the oldest trees worldwide; it is found mainly in the Mediterranean Basin. Genetic pattering studies have shown that cultivated olive trees, i.e., *Olea europaea* L. *var Europaea*, are more similar to oleaster species, providing evidence to support the concept that oleasters are the ancestors of cultivated trees [1]. Both wild and cultivated olive oil have beneficial properties for human health, giving them high economic and nutritional value; however, this has made olive oil one of the most vulnerable agricultural products to fraud and fakery. Wild-type olive oil has higher antioxidant activity, as well as phenolic, tocopherolic and orthodiphenolic contents equal to or higher those in extra virgin cultivated olive oil [2]. Moreover, wild-type olive is a valuable natural resource due to its resistance to certain environmental and climatic conditions and diseases [3]. For the above reasons, its genetic characteristics have to be evaluated, and reliable molecular tools have to be developed for olive oil origin traceability (genetically and geographically) and wild-type olive oil identification. On the other hand, producers need accurate analytical tools for the genetic identification of their wild-type olive-related products to ensure their high added value [4].

Genetic variations between the two plant species have not been extensively explored by the research community. The analytical techniques used so far for the genetic identification of the wild form of olive tree include randomly amplified polymorphic DNA (RAPD), amplified fragment length polymorphisms (AFLPs) and intersimple and simple sequence repeats (ISSRs and SSRs), based on the chloroplastic and mitochondrial plant DNA [1]. Early research compared the genome of *Olea* *europaea* L. *var Europaea* to that of the wild-type olive, derived from many countries and two areas in Italy, using AFLP analysis as designed by Angiolillo et al., 1999, and Baldoni et al., 2006 [5,6]. RAPD analysis was used to distinguish oleasters from *Olea europaea* L. *var Europaea* trees on the Mediterranean islands of Corsica and Sardinia, as well as in Turkey [7,8]. Besnard et al. used RAPD markers and restriction fragment length polymorphism (RFLP) analysis based on mitochondrial and cytoplasmatic DNA to investigate the relationships among olive species and subspecies in the Mediterranean Basin and other countries in Asia and Africa. This research led to the discovery that there was a large degree of diversity among olive cultivated trees, but that they were more or less related to the local oleasters [9,10]. Moreover, ISSR and SSR markers have been utilized by many researchers to investigate the relation and differentiation of cultivated olives from wild-type olives [3,11–16]. Genome size estimation based on double-stranded DNA staining followed by flow cytometric analysis was also used for screening purposes between *Olea europaea var Sylvestris* and *Olea europaea* L. *var Europaea* species [17], while flow cytometry in combination with SSR profiles was used for the taxonomy of four olive subspecies, namely *Olea europaea ssp. cerasiformis*, *Olea europaea ssp. guanchica*, *Olea europaea var Sylvestris* and *Olea europaea* L. *var Europaea* [18].

Moreover, the wild olive has also been used for nonedible purposes in pharmacology and cosmetics to create products with specific valuable characteristics. Researches have also studied the antimicrobial activity of the wild olive against certain human bacterial pathogens [19]. Several plants, including the olive and its wild form, have also been used for the production of various food supplements [20]. Finally, phenolic extracts from wild olive leaves have been investigated for use in foodstuffs, food additives and functional food materials, due to their high antioxidant activity [21,22].

In 2017, the complete genome sequence of *Olea europaea var Sylvestris* was published by Unver et al. [23]. This will be useful, in the future, for the localization of specific genetic variations in the genome of oleasters compared to other olive subspecies.

For the first time, in this work, a single nucleotide polymorphism (SNP)-based method was developed for the detection and identification of the wild form of olive in order to distinguish it from the cultivated olive. Different olive cultivars contain different SNPs in their genome that are responsible for their unique phenotyping characteristics [24,25]. The method was based on an allele-specific, real-time PCR. The proposed method is able to detect wild-type olive DNA at levels as low as 1% in DNA derived from the cultivated olive.

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

#### *2.1. Materials and Instrumentation*

The Vent (exo-) DNA polymerase was purchased by New England Biolabs (Beverly, MA, USA). Deoxynucleoside triphosphates (dNTPs) were obtained from Kapa Biosystems (Wilmington, MA, USA). The fluorescent dye SYBR Green I 104 × concentrated was from Molecular Probes (Eugene, OR, USA). The primers used were from Eurofins Scientific (Brussels, Belgium) and are listed in Table 1. The size of the PCR products was 136 bp. An extra virgin olive oil sample (*Olea europaea* L. *var Europaea*) was purchased from a local market, while a certified wild-type olive oil sample (*Olea europaea var Sylvestris*) was kindly by local producer, Alexandros Karakikes, from the Olea Sylvestris estate (Agrielaio, Volos, Greece) [26].

Real-time PCR was performed using the Mini Opticon Real-Time PCR System from Biorad (Hercules, CA, USA), while the results were analyzed using the Bio–Rad CFX Manager 3.0 software.

**Table 1.** The primers used in the allele-specific, real-time PCR, two species-specific upstream primers and a common downstream primer, along with their melting temperatures (Tm).

\* according to Eurofins Scientific (Brussels, Belgium).

#### *2.2. DNA Isolation Procedure*

DNA was isolated from olive oil samples using the NucleoSpin Tissue kit from Macherey-Nagel (Düren, Germany) according to the manufacturer's instructions. The quantity and purity of the isolated DNA were determined using the Nanodrop UV/VIS Nanophotometer by Implen GmbH (Münich, Germany).

#### *2.3. Design of the Primers*

The primers used for the amplification of *Olea europaea var Sylvestris* (wild-type olive) and *var Europaea* (cultivated olive) were designed using the free online Oligo Analyzer software for primer evaluation (created by Dr. Teemu Kuulasmaa), based on the *Olea europaea var. sylvestris* NADH dehydrogenase subunit F gene, chloroplastic sequence (Accession Number: AY172114) and the *Olea europaea* L. NADH dehydrogenase subunit F (ndhF) gene chloroplastic sequence (Accession Number: DQ673278) [23].

#### *2.4. Allele-Specific, Real-Time PCR*

The allele-specific, real-time PCR reactions were conducted in a final volume of 50 µL and contained 1 × Thermopol Buffer (20 mM Tris-HCl, 10 mM (NH4)2SO4, 10 mM KCl, 0.1% Triton® X-100 at pH 8.8), 0.5 µM of each of the upstream and downstream primers, 0.2 mM of each of the four dNTPs, 0.5 mM MgCl2, 2 × SYBR Green I, one unit of Vent (exo-) DNA polymerase and 150 ng of isolated DNA. The reaction conditions involved a 95 ◦C incubation step for three min, followed by 45 cycles at 95 ◦C for 30 s, 62 ◦C for 30 s, 72 ◦C for 30 s and a final extension step at 72 ◦C for 10 min.

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

A new analytical method was developed for the detection and identification of *Olea europaea var Sylvestris* that refers to the wild form of the olive tree. The method was based on the detection of a specific Single Nucleotide Polymorphism (SNP) that is different in the genome of the wild olive plant. The method involves the following steps: (i) DNA isolation from olive oil samples and (ii) allele-specific, real-time PCR using an upstream primer specific to *Olea europaea var Sylvestris* or *var Europaea* species and a common downstream primer. The species-specific primers have the same 22-base sequence but differ only at the base at the 30 end that contains the SNP of interest. The DNA sequences were amplified using a DNA polymerase that lacked the 30 to 50 exonuclease activity, so only the primer that was perfectly complementary to the DNA target was extended by the enzyme. The amplicons were finally detected using the DNA intercalating fluorescent dye SYBR Green I. The principle of the proposed method is illustrated in Figure 1. SYBR Green I was chosen here instead of Taqman probes in order to develop a new analytical method that could be easily transferred, with few modifications, for the detection of other SNPs that will be found in the wild olive genome in the future.

*Foods* **2020**, *9*, x FOR PEER REVIEW 4 of 9

**Figure 1. (Upper panel)** Schematic illustration of the principle of the method that includes DNA isolation and purification from olive oil samples using spin cleanup columns, including the following steps: cell lysis of an olive oil sample, capture of DNA to the cleanup columns and elution of the DNA from the columns. **(Lower panel)** The allele-specific, real-time PCR. Two allele-specific upstream primers that contain the SNP of interest at their 3΄ ends and one common downstream primer were used in the amplification reaction. Only the perfectly complementary upstream primer to the target was extended by the DNA polymerase, while the amplicons were detected by the DNA intercalating dye, SYBR Green I. *3.1. DNA Isolation*  **Figure 1.** (**Upper panel**) Schematic illustration of the principle of the method that includes DNA isolation and purification from olive oil samples using spin cleanup columns, including the following steps: cell lysis of an olive oil sample, capture of DNA to the cleanup columns and elution of the DNA from the columns. (**Lower panel**) The allele-specific, real-time PCR. Two allele-specific upstream primers that contain the SNP of interest at their 30 ends and one common downstream primer were used in the amplification reaction. Only the perfectly complementary upstream primer to the target was extended by the DNA polymerase, while the amplicons were detected by the DNA intercalating dye, SYBR Green I.

#### First, DNA was isolated from olive oil samples and its concentration was determined using a *3.1. DNA Isolation*

UV/VIS nanophotometer. It was found that the isolation procedure did not result in a constant DNA amount for all samples, with the DNA concentrations ranging from 8.4 to 142 ng/μL. To avoid fluctuation in the PCR yield due to different initial DNA concentrations, we decided to use the same amount (ng) of isolated DNA for all samples into the real-time PCR mixture. After amplification, the amplicons had a size of 136 bp. The quality of the isolated DNA was also determined by UV measurements; the ratios A260/A280 were from 1,174 to 1,739. DNA was considered to be of high quality when the ratio A260/A280 was above 1.8. *3.2. Optimization of the PCR Conditions*  The real-time PCR conditions were initially optimized. The parameters studied were the amount First, DNA was isolated from olive oil samples and its concentration was determined using a UV/VIS nanophotometer. It was found that the isolation procedure did not result in a constant DNA amount for all samples, with the DNA concentrations ranging from 8.4 to 142 ng/µL. To avoid fluctuation in the PCR yield due to different initial DNA concentrations, we decided to use the same amount (ng) of isolated DNA for all samples into the real-time PCR mixture. After amplification, the amplicons had a size of 136 bp. The quality of the isolated DNA was also determined by UV measurements; the ratios A260/A280 were from 1174 to 1739. DNA was considered to be of high quality when the ratio A260/A280 was above 1.8.

#### of the isolated DNA, the concentration of the primers, the number of PCR cycles and the temperature of the annealing step of the reaction. At low DNA and primer concentrations, low temperature (55- *3.2. Optimization of the PCR Conditions*

60 °C) and number of cycles < 45, the PCR was not sufficiently efficient. The yield of the reaction also decreased when a high amount of initial DNA target was used. This may be attributed to the fact that the DNA isolated from olive samples has reduced quality, as it contains high amounts of PCR inhibitors that may inhibit the activity of the DNA polymerase [27]. We also observed that the highest reaction yield and specificity were obtained at an annealing temperature of 62 °C. The real-time PCR conditions were initially optimized. The parameters studied were the amount of the isolated DNA, the concentration of the primers, the number of PCR cycles and the temperature of the annealing step of the reaction. At low DNA and primer concentrations, low temperature (55–60 ◦C) and number of cycles < 45, the PCR was not sufficiently efficient. The yield of the reaction also decreased when a high amount of initial DNA target was used. This may be attributed to the fact that the DNA isolated from olive samples has reduced quality, as it contains high amounts of PCR inhibitors that may inhibit the activity of the DNA polymerase [27]. We also observed that the highest reaction yield and specificity were obtained at an annealing temperature of 62 ◦C.

#### *3.3. Specificity of the Allele-Specific Primers Foods* **2020**, *9*, x FOR PEER REVIEW 5 of 9

The specificity of the two species-dependent upstream primers was then studied as follows: both DNA targets, *Olea europaea var Sylvestris* (wild-type olive) and *var Europaea* (cultivated olive) were subjected to two separate amplification reactions using either the upstream primer specific to the wild-type olive or the cultivated olive-specific upstream primer. As shown, in Figure 2, each primer amplified only its fully complementary DNA sequence, proving the superior specificity of the primers. To ensure that the fluorescence signals were attributed only to the specific amplicons, a melting curve analysis was also performed after each amplification reaction. The melting curve analysis revealed only one peak for each PCR product, the melting temperature (Tm) of which was 77 ◦C for *Olea europaea var Sylvestris* (wild-type olive) and 78 ◦C for *var Europaea* (cultivated olive), allowing us to distinguish between the two allele-specific DNA sequences. *3.3. Specificity of the Allele-Specific Primers*  The specificity of the two species-dependent upstream primers was then studied as follows: both DNA targets, *Olea europaea var Sylvestris* (wild-type olive) and *var Europaea* (cultivated olive) were subjected to two separate amplification reactions using either the upstream primer specific to the wild-type olive or the cultivated olive-specific upstream primer. As shown, in Figure 2, each primer amplified only its fully complementary DNA sequence, proving the superior specificity of the primers. To ensure that the fluorescence signals were attributed only to the specific amplicons, a melting curve analysis was also performed after each amplification reaction. The melting curve analysis revealed only one peak for each PCR product, the melting temperature (Tm) of which was 77 °C for *Olea europaea var Sylvestris* (wild-type olive) and 78 °C for *var Europaea* (cultivated olive),

allowing us to distinguish between the two allele-specific DNA sequences.

**Figure 2.** The real-time PCR curves, along with the corresponding melting curve analysis, obtained during the specificity study of the two-allele specific upstream primers with both DNA targets: *Olea europaea var Sylvestris* (wild-type of olive) **(a)** and *Olea europaea* L*. var Europaea* (cultivated olive) **(b)**. Each specific primer strictly amplifies the fully complementary DNA sequence. Tm: melting temperature, RFU: Relative Fluorescence Units. **Figure 2.** The real-time PCR curves, along with the corresponding melting curve analysis, obtained during the specificity study of the two-allele specific upstream primers with both DNA targets: *Olea europaea var Sylvestris* (wild-type of olive) (**a**) and *Olea europaea* L. *var Europaea* (cultivated olive) (**b**). Each specific primer strictly amplifies the fully complementary DNA sequence. Tm: melting temperature, RFU: Relative Fluorescence Units.

#### *3.4. Detectability of the Method in Binary DNA Mixtures*

Subsequently, the detectability of the method in olive DNA binary mixtures was evaluated. DNA mixtures that contained different proportions (1–50%) of DNA from *Olea europaea var Sylvestris* in DNA from *var Europaea* were prepared. An amount of 150 ng of each DNA mixture was then subjected to two separate allele-specific, real-time PCR reactions using each of the species-specific upstream primers along with the common downstream primer, respectively. A high amount of total DNA was used in the PCR in order to detect the low amount of wild olive DNA in the mixtures, e.g., for 150 ng of total DNA in the 1% mixture, only the 1.5 ng was the wild olive DNA. The results are presented in Figure 3. We were able to detect as little as 1% of DNA specific to *Olea europaea var Sylvestris* in the presence of DNA from *Olea europaea* L. *var Europaea*. The allelic ratios of the analyzed SNP for the above DNA mixtures were also calculated based on the fluorescence value at the 45th cycle of the reaction, and are

presented in the same Figure. The allelic ratios for all DNA mixtures were close to the value of 0.5, as expected for a heterozygote sample. the analyzed SNP for the above DNA mixtures were also calculated based on the fluorescence value at the 45th cycle of the reaction, and are presented in the same Figure. The allelic ratios for all DNA mixtures were close to the value of 0.5, as expected for a heterozygote sample.

*europaea var Sylvestris* in the presence of DNA from *Olea europaea* L. *var Europaea*. The allelic ratios of

*Foods* **2020**, *9*, x FOR PEER REVIEW 6 of 9

Subsequently, the detectability of the method in olive DNA binary mixtures was evaluated. DNA mixtures that contained different proportions (1%–50%) of DNA from *Olea europaea var Sylvestris* in DNA from *var Europaea* were prepared. An amount of 150 ng of each DNA mixture was then subjected to two separate allele-specific, real-time PCR reactions using each of the speciesspecific upstream primers along with the common downstream primer, respectively. A high amount of total DNA was used in the PCR in order to detect the low amount of wild olive DNA in the mixtures, e.g., for 150 ng of total DNA in the 1% mixture, only the 1.5 ng was the wild olive DNA.

*3.4 Detectability of the Method in Binary DNA Mixtures*

**Figure 3. (a)** The real-time PCR curves obtained for different % DNA content (0%–50%) of *Olea europaea var Sylvestris* (wild-type olive) DNA in binary mixtures with *Olea europaea* L*. var Europaea* (cultivated olive) DNA. **(b)** The allelic ratios of the binary DNA mixtures calculated as the ratio of the fluorescence intensity obtained with the upstream primer specific to *Olea europaea var Sylvestris* target versus the sum of the fluorescence intensity obtained by both allele-specific primers for *Olea europaea var Sylvestris* and *var Europaea* targets. All allelic ratios were close to the value of 0.5, which corresponds to a heterozygote sample. RFU: Relative Fluorescence Units. **Figure 3.** (**a**) The real-time PCR curves obtained for different % DNA content (0–50%) of *Olea europaea var Sylvestris* (wild-type olive) DNA in binary mixtures with *Olea europaea* L. *var Europaea* (cultivated olive) DNA. (**b**) The allelic ratios of the binary DNA mixtures calculated as the ratio of the fluorescence intensity obtained with the upstream primer specific to *Olea europaea var Sylvestris* target versus the sum of the fluorescence intensity obtained by both allele-specific primers for *Olea europaea var Sylvestris* and *var Europaea* targets. All allelic ratios were close to the value of 0.5, which corresponds to a heterozygote sample. RFU: Relative Fluorescence Units.

#### *3.5 Reproducibility of the Method 3.5. Reproducibility of the Method*

Finally, the reproducibility of the method was determined. Two different proportions, 1% and 10%, of the above DNA mixtures, were subjected, in triplicate, to real-time PCR. The % coefficients of variation (CV) were calculated based on the obtained Cq values for all samples. The CV for the 1%-content was 10.5% and for the 10%-content was 7.5%, demonstrating the reproducibility of the method. Finally, the reproducibility of the method was determined. Two different proportions, 1% and 10%, of the above DNA mixtures, were subjected, in triplicate, to real-time PCR. The % coefficients of variation (CV) were calculated based on the obtained Cq values for all samples. The CV for the 1%-content was 10.5% and for the 10%-content was 7.5%, demonstrating the reproducibility of the method.

#### **4. Conclusions 4. Conclusions**

A new allele-specific, real-time PCR-based analytical method was developed for the detection and identification of wild-type olive oil (*Olea europaea var Sylvestris*), compared to cultivated olive oil (*Olea europaea L. var Europaea*). The discrimination of the two similar plant species was based on genotyping a single nucleotide polymorphism (SNP) that is differently present in the genome of the two plant species. The detection of this SNP was carried out by an allele-specific, real-time PCR that was performed using two different species-specific upstream primers that contained the analyzed SNP and a common downstream primer. Each specific primer amplified only its fully complementary DNA sequence, leading to species identification. The detection of the amplicons was accomplished A new allele-specific, real-time PCR-based analytical method was developed for the detection and identification of wild-type olive oil (*Olea europaea var Sylvestris*), compared to cultivated olive oil (*Olea europaea* L. *var Europaea*). The discrimination of the two similar plant species was based on genotyping a single nucleotide polymorphism (SNP) that is differently present in the genome of the two plant species. The detection of this SNP was carried out by an allele-specific, real-time PCR that was performed using two different species-specific upstream primers that contained the analyzed SNP and a common downstream primer. Each specific primer amplified only its fully complementary DNA sequence, leading to species identification. The detection of the amplicons was accomplished using the DNA intercalating dye, SYBR Green I. With the proposed method, we were able to sucessfully distinguish between the two plant species in olive oil samples. Also, as little as 1% wild-type olive species was detected in binary DNA mixtures of the two analyzed plant species. In conclusion, the method is easy, rapid, has good detectability, is reproducible and can easily distinguish between species. The proposed method also contributes to the ability to add the higher financial value to wild-type olive-based products. In the future, the determination of different SNPs in the wild-type olive genome compared to all the known cultivated olive trees could lead to more accurate discrimination of the wild-type olive among other olive-based subspecies. The proposed method could also be applied, with some modifications, for the detection of wild olive-based ingredients in food supplements and cosmetic products. The global increase in food supplements has led to the mislabeling of these products and fraudulent practices. In both cases, the purity of the extracted DNA is more important than the

PCR yield itself, because several food additives and other ingredients may be present in the extracts, inhibiting the PCR amplification. Also, the amount of the extracted DNA may be extremely low. Thus, the DNA isolation protocols have to be properly justified to remove these inhibitors and increase the DNA recovery and the PCR yield. In some studies, however, the inability to extract DNA from some food supplements has been reported. Finally, in some products, DNA degradation may also occur due to thermal or chemical treatment, but the use of short-length amplicons can overcome this issue [28–31].

**Author Contributions:** Conceptualization, D.P.K.; methodology, C.I.K. and D.P.K.; software, C.I.K.; investigation, C.I.K. and D.P.K.; resources, C.I.K. and D.P.K.; data curation, C.I.K. and D.P.K.; writing—original draft preparation, C.I.K. and D.P.K.; writing—review and editing, D.P.K.; supervision, D.P.K.; project administration, D.P.K.; funding acquisition, D.P.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the project "Research Infrastructure on Food Bioprocessing Development and Innovation Exploitation—Food Innovation RI" (MIS 5027222), which is implemented under the Action "Reinforcement of the Research and Innovation Infrastructure", funded by the Operational Programme "Competitiveness, Entrepreneurship and Innovation" (NSRF 2014–2020) and co-financed by Greece and the European Union (European Regional Development Fund).

**Acknowledgments:** We acknowledge Alexandros Karakikes (Olea Sylvestris estate) for kindly providing us with the wild-type olive oil sample.

**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*

## **Using Chemometric Analyses for Tracing the Regional Origin of Multifloral Honeys of Montenegro**

**Vesna Vukašinovi´c-Peši´c 1,\*, Nada Blagojevi´c 1 , Snežana Brašanac-Vukanovi´c 1 , Ana Savi´c 2 and Vladimir Peši´c 3**

Received: 19 January 2020; Accepted: 14 February 2020; Published: 18 February 2020

**Abstract:** This is the first study of mineral content and basic physicochemical parameters of honeys of Montenegro. We examined honey samples from eight different micro-regions of Montenegro, and the results confirm that, with the exception of cadmium in samples from two regions exposed to industrial pollution, none of the 12 elements analyzed exceeded the maximum allowable level. The samples from areas exposed to industrial pollution were clearly distinguished from samples from other regions of Montenegro in the detectable contents of Pb, Cd, and Sr. This study showed that chemometric techniques might enhance the classification of Montenegrin honeys according to their micro-regional origin using the mineral content. Linear discriminant analysis revealed that the classification rate was 79.2% using the cross-validation method.

**Keywords:** honey; regional origin; chemometric analysis; mineral content; Montenegro

### **1. Introduction**

Honey is a complex natural product, whose characteristics depend on the flower nectar from which it is obtained, but also on other factors such as geographical origin, bee species, season, type of processing and storage [1]. It is known that pollution and a number of different pollutants present in its foraging areas have an impact on honeybees [2] but also on nectar-providing plant species. Therefore, it is necessary to assure geographical traceability and determine the botanical origin of the foraging area of the beehive.

As stated by Karabagias and Karabournioti [3], the authentication of honey is gaining in importance and includes a number of contending parties from producers and sellers to consumers and control labs. A number of papers have shown that specific physicochemical parameters and mineral contents in combination with chemometric analyses can be a useful tool in discovering botanical and/or geographical origin of honeys that may enter the market [1,3,4].

Tracing the geographical origin of honeys can provide important information about the potential contamination of the area from which the honey production material comes. Therefore, ensuring high standards in terms of product safety leads to the need to examine the contents of essential and toxic elements in honey. Due to its bioaccumulation ability, honey can be used as an indicator of metal pollution, especially of toxic pollutants such as Pb, Cd, and As [5–7].

Due to its geographical position, climate conditions and richness of the nectar-providing plants, Montenegro provides favorable natural conditions for more intensive development of beekeeping. According to the data for 2011, the population in Montenegro was 625,266, while the honey production

for that year was 394 t, and the average annual consumption 1.2 kg per person, meaning that a large part of honey consumption in Montenegro is imported [8]. Data for the last few years show an increase in honey production (627 t for 2016) but also an increase in the average annual consumption of honey per person (2.76 kg) [9].

The majority of honey on the market in Montenegro are multifloral (derived from a large number of nectar-providing plant species in the honeybees foraging area). Most of these honey types are recognizable by their local or regional origin (e.g., Katunski med (= honey), Pivski med, Piperski med.). It is worth mentioning that Montenegro and its regions are known to harbor a high number of regional floral endemics [10] that likely affect the composition and properties of honey.

There is a lack of information on the mineral content and basic physicochemical parameters of honey from the territory of Montenegro. Moreover, there is no continuity in monitoring the quality of honey, especially in areas that are exposed to the effects of potential pollution sources. Due to the high consumption of local honey in the diet, the need and obligation for its systematic characterization are highly required.

This study is aimed to investigate the mineral content and the basic physicochemical parameters of honeys from different micro-regions of Montenegro. We evaluated the usefulness of chemometric analyses for the classification of honeys according to its regional origin.

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

Twenty-four honey samples as indicated in Figure 1 were collected from eight micro-regions of Montenegro, i.e., (1) Piva, (2) Zbljevo, (3) Potrlica, (4) Mijakovi´ci, (5) Piperi, (6) Martini´ci, (7) Katunska, and (8) Zeta. The Piva, Zbljevo, Potrlica and Mijakovi´ci micro-regions are situated in the continental part of Montenegro (Alpine biogeographical region, see Figure 1) while the four other micro-regions are situated in the sub-Mediterranean part of the country belonging to the Mediterranean biogeographical region [10]. The climate in the latter region is mainly Mediterranean-Adriatic with relatively dry and warm summers (the average air temperature of the warmest month > 20 ◦C), but humid and mild winters (the average air temperatures varies from 6 to 9 ◦C), while the Alpine region has a "continental" type climate, with relatively cool and humid summers and long and harsh winters [10].

Samples were taken from individual beekeepers during the harvesting season 2015. All samples were multifloral as confirmed by the suppliers. The samples were stored in glass flasks at room temperature before analysis. Physicochemical parameters (pH, electrical conductivity (EC), free acidity (FA) and moisture) were analyzed using the Harmonized Methods of the International Honey Commission [11].

The mineral composition of honey was analyzed by inductively coupled plasma-optical emission spectrometry (ICP-OES). About 1 g of each honey sample was digested with 14 mL 65% HNO3 and 2 mL 35% H2O2 on a hot plate to near dryness. The sample containing a volumetric flask was cooled at room temperature before the addition of deionized water to the mark on the flask. All samples were prepared in triplicate and their average value was assessed.

The concentration of twelve elements (Pb, Cd, Cu, Zn, Fe, Cr, Sr, Ba, Ca, Na, K, Mg) were determined by ICP-OES according to the iCAP 6000 spectrometer method.

All statistical analyses were performed using SPSS 17.0 (SPSS Statistics for Windows, Version 17.0. SPSS Inc., Chicago, IL, USA). Data were expressed as mean ± standard deviation. A Kolmogorov–Smirnov test showed that all analyzed physicochemical parameters were normally distributed, while the content of Pb, Cd, Sr and Ba in some regions exhibited significant differences from the normal distribution. The one-way analysis of variance (ANOVA) was performed on physicochemical parameters in order to determine if there any significant differences between studied micro-regions at the confidence level 0.05. The Kruskal–Wallis test was used to investigate whether the mineral contents varied significantly between the investigated micro-regions. The relationship between the mineral content and physicochemical parameters were analyzed using the Spearman's correlation analysis. For checking similarities between samples of honey of different geographical

origin we used two chemometric analyses: PCA and LDA. Statistical methods based on principal component analysis (PCA) and linear discriminant analysis (LDA) have been used. The LDA was performed using R.3.5.3, while the PCA was made by using MVSP version 3.21. *Foods* **2020**, *9*, x FOR PEER REVIEW 3 of 9

**Figure 1.** Map of Montenegro with marked locations of honey sampling from eight micro‐regions (in parentheses are given sampling location numbers): Piva (1–3), Zbljevo (4–6), Potrlica (7–9), Mijakovići (10–12), Piperi (13–15), Martinići (16–18), Katunska (19–21), and Zeta (22–24). **Figure 1.** Map of Montenegro with marked locations of honey sampling from eight micro-regions (in parentheses are given sampling location numbers): Piva (1–3), Zbljevo (4–6), Potrlica (7–9), Mijakovi´ci (10–12), Piperi (13–15), Martini´ci (16–18), Katunska (19–21), and Zeta (22–24).

#### All statistical analyses were performed using SPSS 17.0 (SPSS Statistics for Windows, Version **3. Results**

**3. Results**

Zbljevo.

17.0. SPSS Inc., Chicago, IL, USA). Data were expressed as mean ± standard deviation. A Kolmogorov–Smirnov test showed that all analyzed physicochemical parameters were normally distributed, while the content of Pb, Cd, Sr and Ba in some regions exhibited significant differences from the normal distribution. The one‐way analysis of variance (ANOVA) was performed on physicochemical parameters in order to determine if there any significant differences between studied micro‐regions at the confidence level 0.05. The Kruskal–Wallis test was used to investigate whether the mineral contents varied significantly between the investigated micro‐regions. The relationship between the mineral content and physicochemical parameters were analyzed using the The mineral content of honey samples from different geographical areas of Montenegro is presented in Table 1. The value presented for each element is the average concentration observed. A significant difference has been observed in the concentrations of Pb, Cd and Sr (*p* = 0.002) between studied micro-regions. In most analyzed samples the concentrations of above-listed elements were below the limit of detection except in the samples from Potrlica, Zbljevo and Mijakovi´ci. The highest Cd concentration was observed in samples from Potrlica (0.08 ± 0.01 mg/kg). The highest concentration of Pb (0.21 ± 0.06 mg/kg) and Sr (0.12 ± 0.00 mg/kg) were recorded in samples from Zbljevo.

Spearman's correlation analysis. For checking similarities between samples of honey of different geographical origin we used two chemometric analyses: PCA and LDA. Statistical methods based on principal component analysis (PCA) and linear discriminant analysis (LDA) have been used. The

LDA was performed using R.3.5.3, while the PCA was made by using MVSP version 3.21.

The mineral content of honey samples from different geographical areas of Montenegro is presented in Table 1. The value presented for each element is the average concentration observed. A significant difference has been observed in the concentrations of Pb, Cd and Sr (*p* = 0.002) between studied micro‐regions. In most analyzed samples the concentrations of above‐listed elements were

*Foods* **2020**, *9*, 210

**Table 1.** Mineral content and physicochemical parameters of honey samples from studied micro-regions of Montenegro.

72

The concentrations of examined physicochemical parameters of honey samples are given in Table 1. The moisture level had similar values across studied regions and ranged from 14.92 ± 0.78% (Piperi) to 16.22 ± 0.36% (Zeta). The pH of studied honey samples varies between 3.87 and 4.49 and was lowest in samples from the Piva region (3.87 ± 0.36) and highest in honey samples from Mijakovi´ci (4.49 ± 0.14). A significant difference has been observed in pH according to the honey regional origin (*p* = 0.048). The electrical conductivity varied from 0.39 to 0.93 mS/cm and was lowest in samples from Katunska (0.39 ± 0.08 mS/cm) and highest in honey samples from Zeta (0.93 ± 0.15 mS/cm). A significant difference has been observed in electrical conductivity according to the honey geographical origin (*p* = 0.013). Free acidity varied from 25.00 to 41.67 meq/kg and was lowest in samples from Katunska (25.00 ± 7.21 meq/kg) and highest in honey samples from the Piva region (41.67 ± 12.10 meq/kg).

Correlation analysis revealed significant correlation between contents of K (R = 0.800, significance < 0.001) and Mg (R = 0.758, significance < 0.001) from one side and pH from the other one (Table 2).

**Table 2.** Results of the correlation analysis between the mineral content and physicochemical parameters of analyzed honey samples from Montenegro.

\*\* significance < 0.01.

The first principal component explains 42.53% of the total variability and is mostly determined by Cd (R = 0.416), Pb (R = 0.398), and Cu (R = 0.352). The PC2 explains 17.97% and is mostly determined by Mg (R = −0.577), K (R = −0.54), and Na (R = −0.425). Mutual projections of factor scores and their loadings for the first two PCs are presented in Figure 2. As can be seen from the projection plot the separation of the analyzed honey samples is much clearer along the *X*-axis. On the one side, there are localities Piva, Piperi, Katunska, Zeta and Martini´ci in whose honey samples Cd, Pb, and Sr were not detected. On the other side, there are Mijakovi´ci, Potrlica and especially Zbljevo, whose honey samples concentrations of Cd, Pb and Sr were detected.

73

*Foods* **2020**, *9*, x FOR PEER REVIEW 6 of 9

**Figure 2.** Principal component analysis (PCA) of the mineral content scores of analyzed Montenegrin honey samples. **Figure 2.** Principal component analysis (PCA) of the mineral content scores of analyzed Montenegrin honey samples.

LDA performed on the geographical origin revealed that the cross‐validation classification was correct for 79.17% of samples (Table 3). The smallest percent of good classification was achieved in the case of honey samples from Katunska, while the highest in the case of honeys from Piperi, Zeta, Mijakovići, and Zbljevo. LDA performed on the geographical origin revealed that the cross-validation classification was correct for 79.17% of samples (Table 3). The smallest percent of good classification was achieved in the case of honey samples from Katunska, while the highest in the case of honeys from Piperi, Zeta, Mijakovi´ci, and Zbljevo.

**Table 3.** Classification of honey according to their regional origin using the linear discriminate **Table 3.** Classification of honey according to their regional origin using the linear discriminate analysis.

Zbljevo 0.00% 0.00% 33.33% 0.00% 0.00% 0.00% 33.33% 100.00% The bidimensional plot (Figure 3) of the first two functions show four distinct clusters, three of them corresponding to Mijakovići, Potrlica, and Zbljevo regions, while all other regions were clustered together. The first discriminant function explains 94.2% of the total variance and it is The bidimensional plot (Figure 3) of the first two functions show four distinct clusters, three of them corresponding to Mijakovi´ci, Potrlica, and Zbljevo regions, while all other regions were clustered together. The first discriminant function explains 94.2% of the total variance and it is dominated by Cd content (R = 0.95). The second discriminant function explains 4.4% of the total variance and is dominated by Pb (R = −0.55) and Sr (R = −0.59) contents.

dominated by Cd content (R = 0.95). The second discriminant function explains 4.4% of the total

variance and is dominated by Pb (R = −0.55) and Sr (R = −0.59) contents.

Potrlica 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 66.67% 0.00%

*Foods* **2020**, *9*, x FOR PEER REVIEW 7 of 9

**Figure 3.** Linear discriminant score plot of analyzed honey samples. **Figure 3.** Linear discriminant score plot of analyzed honey samples.

#### **4. Discussion 4. Discussion**

The values of the mineral contents have been compared with those established by the EU regulations [12]. With the exception of the concentrations of cadmium in samples from Zbljevo and Potrlica, none of the 12 elements analyzed exceeded the maximum allowable level established by the EU regulations. Our results revealed significant differences in the concentrations of Pb, Cd, and Sr between the studied geographical areas of Montenegro. The latter elements were detected only in the samples from Mijakovići, Potrlica and Zbljevo regions, which are likely under the influence of the Pljevlja Thermal Power Plant (Zbljevo, and in less extent Mijakovići) and the Pljevlja coalmine (Potrlica). The values of the mineral contents have been compared with those established by the EU regulations [12]. With the exception of the concentrations of cadmium in samples from Zbljevo and Potrlica, none of the 12 elements analyzed exceeded the maximum allowable level established by the EU regulations. Our results revealed significant differences in the concentrations of Pb, Cd, and Sr between the studied geographical areas of Montenegro. The latter elements were detected only in the samples from Mijakovi´ci, Potrlica and Zbljevo regions, which are likely under the influence of the Pljevlja Thermal Power Plant (Zbljevo, and in less extent Mijakovi´ci) and the Pljevlja coalmine (Potrlica).

The Pb content in the analyzed honey samples varied from 80–210 μg/kg. These values were lower in comparison with the honey from Serbia (290 μg/kg [13]) and Italy (289 μg/kg [14]) but higher in comparison with those from Croatia (5.43–11.3 μg/kg [15]) and Bosnia and Hercegovina (13.4 μg/kg [16]). All these values are below the maximum allowable level established by the EU regulations (0.5 mg/kg) [12]. The Pb content in the analyzed honey samples varied from 80–210 µg/kg. These values were lower in comparison with the honey from Serbia (290 µg/kg [13]) and Italy (289 µg/kg [14]) but higher in comparison with those from Croatia (5.43–11.3 µg/kg [15]) and Bosnia and Hercegovina (13.4 µg/kg [16]). All these values are below the maximum allowable level established by the EU regulations (0.5 mg/kg) [12].

The values reported for Cd in this study (20–80 μg/kg) were higher in comparison with the honey from Croatia (0.69–12.8 μg/kg [15]), Bosnia and Hercegovina (0.013–22.9 μg/kg [16]), Romania (0.5–11.60 μg/kg [17]), Italy (8–18 μg/kg [14]), Spain (0.7–50 μg/kg [17]) and Serbia (0.59–30 μg/kg [1,13]). The values of Cd content in the samples from Potrlica and Zbljevo exceeds the maximum allowable level established by the EU legislation (0.05 mg/kg) [12]. As the main sources of Cd are recognized as the presence in sewage sludge and smelting from the nearby Pljevlja Thermal Power Plant (Zbljevo), or mining from the Pljevlja coalmine (Potrlica). The values reported for Cd in this study (20–80 µg/kg) were higher in comparison with the honey from Croatia (0.69–12.8 µg/kg [15]), Bosnia and Hercegovina (0.013–22.9 µg/kg [16]), Romania (0.5–11.60 µg/kg [17]), Italy (8–18 µg/kg [14]), Spain (0.7–50 µg/kg [17]) and Serbia (0.59–30 µg/kg [1,13]). The values of Cd content in the samples from Potrlica and Zbljevo exceeds the maximum allowable level established by the EU legislation (0.05 mg/kg) [12]. As the main sources of Cd are recognized as the presence in sewage sludge and smelting from the nearby Pljevlja Thermal Power Plant (Zbljevo), or mining from the Pljevlja coalmine (Potrlica).

The Sr content in honey samples from our study varied from 0.07–0.12 μg/kg and was in the same range as those from Serbia (0.09–0.19 μg/kg [13]). The Sr content in honey samples from our study varied from 0.07–0.12 µg/kg and was in the same range as those from Serbia (0.09–0.19 µg/kg [13]).

The most abundant element in honey samples was kalium, followed by Ca, Mg, Na, and Fe. In our study, we found that the content of kalium and magnesium correlated with pH. The average

The most abundant element in honey samples was kalium, followed by Ca, Mg, Na, and Fe. In our study, we found that the content of kalium and magnesium correlated with pH. The average levels of K content ranged from 713–2589.33 mg/kg and was in the same range with those from Croatia (304.7–2824.4 mg/kg [15]), but lower in comparison with the maximum values established for honeys from Bosnia and Herzegovina (14.81–4895.73 µg/kg [16]). On the other hand, the range of concentrations of kalium in the honey from Serbia (400–1755 mg/kg [1,13]) and Slovenia (1090–1220 mg/kg [18]) were lower. The Mg content in honey samples from our study varied from 29.52 to 76.33 mg/kg. In neighboring countries the Mg content in honey varied in a similar range, as: 28.83 to 101.50 mg/kg [13] in Serbia, 2.18 to 166.04 mg/kg [16] in Bosnia and Herzegovina and from 8.02 to 59.1 mg/kg [16,19] in Croatia.

In our study, we used two chemometric analyses, PCA and LDA, respectively to test similarities between honey samples of honey of different geographical origins. Both applied methods separated the regions exposed to industrial pollution (Mijakovi´ci, Potrlica, and Zbljevo) which are characterized by detectable content of Cd, Pb and Sr in their honey samples.

Using LDA it's possible to evaluate the capacity to correctly predict the group to which the unknown samples belong. In our study LDA analysis performed on the geographical origin revealed that the cross-validation classification was correct for 79.17% of the samples. The obtained values are in the range for those from Serbia (Zlatibor: 94.73%, Vojvodina: 70.58% [1]). On the other hand, our value was greater than those reported in the case of Romania honeys where only 21.2% were correctly classified according to their geographical origin [4].

The smallest percentage of good classification was achieved in the case of honeys from Katunska. Of the three samples from the latter region, only one was correctly classified, while the other two being misclassified as Piperi and Zbljevo, respectively. It is known that large numbers of beekeepers (especially from Katunska) in a part of the year (most often in summertime) move their bee colonies to geographically distant areas. On the other hand, the highest percentage of good classification was achieved in the case of honeys from Piperi, Zeta, Mijakovi´ci, and Zbljevo. One cause may be that most of these sites (i.e., Zeta, Mijakovi´ci, and Zbljevo) are more exposed to industrial pollution, resulting in increased concentration of heavy metals (Pb, Cd, and Sr showing significant difference (*p* < 0.05) between studied regions) in their honeys, which, in turn, increase the success rate of the classification of honey according to their geographical origin.

**Author Contributions:** Investigation, V.V.-P., S.B.-V. and N.B.; statistical analysis, A.S.; writing—original draft preparation, V.P.; writing—review and editing, V.V.-P. and V.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**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*

## **Geographical Origin Discrimination of Monofloral Honeys by Direct Analysis in Real Time Ionization-High Resolution Mass Spectrometry (DART-HRMS)**

#### **Vincenzo Lippolis 1,\* , Elisabetta De Angelis 1 , Giuseppina Maria Fiorino 1 , Annalisa Di Gioia 1 , Marco Arlorio 2 , Antonio Francesco Logrieco 1 and Linda Monaci 1**

### Received: 27 July 2020; Accepted: 28 August 2020; Published: 1 September 2020

**Abstract:** An untargeted method using direct analysis in real time and high resolution mass spectrometry (DART-HRMS) combined to multivariate statistical analysis was developed for the discrimination of two monofloral (chestnut and acacia) honeys for their geographical origins—i.e., Italy and Portugal for chestnut honey and Italy and China for acacia honey. Principal Component Analysis, used as an unsupervised approach, showed samples of clusterization for chestnut honey samples, while overlapping regions were observed for acacia honeys. Three supervised statistical approaches, such as Principal Components—Linear Discriminant Analysis, Partial Least Squares—Discriminant Analysis and k-nearest neighbors, were tested on the dataset gathered and relevant performances were compared. All tested statistical approaches provided comparable prediction abilities in cross-validation and external validation with mean values falling between 89.2–98.4% for chestnut and between 85.8–95.0% for acacia honey. The results obtained herein indicate the feasibility of the DART-HRMS approach in combination with chemometrics for the rapid authentication of honey's geographical origin.

**Keywords:** monofloral honey; direct analysis in real time (DART); high resolution mass spectrometry (HRMS); geographical origin; chemometrics

#### **1. Introduction**

Honey is a complex and high-quality natural product containing a wide range of nutritional and therapeutic properties but with a limited production and high commercial prices. Honey is defined by European Union legislation as the natural sweet substance produced by bees of *Apis mellifera* species from nectar or sugary secretions of plants, as well as from excretions of plant-sucking insects on the living parts of plants [1]. Both the European Union and Codex Alimentarius laws establish that the geographical origin, in terms of country of production, must be indicated on the label, also supplemented by specific reference to the floral or vegetable origin. Moreover, in the case of blends of honey, their origin should be declared as a "blend of EC honeys", "blend of non-EC honeys" or "blend of EC and non-EC honeys" [1,2].

Geographical and botanical origins of honey account for the peculiar chemical composition and organoleptic characteristics of the final product [3,4]. Monofloral honeys, mostly deriving from a single plant species (at least 45% of pollen grains), may considerably differ in their sensory properties with highly prominent flavor and aroma. Acacia (*Robinia pseudoacacia*) honey is one of the most consumed monofloral honeys in Europe, being appreciated for its permanently liquid state, light color, floral aroma and sweet and delicate taste [5]. Similarly, chestnut (*Castanea sativa*) honey is considered one of the most delicious and high-quality honeys, being a very good source for nectar and pollen [6,7]. For these reasons, monofloral honeys, and in particular those derived from acacia and chestnut, have recently gained consumer preferences, with an increased demand and commercial value [4]. Due to their increased commercial value, monofloral honeys are highly susceptible to fraudulent practices through mislabeling and mixing with cheaper and lower-quality honeys or with various sugar syrups.

Honey is produced in different areas of the world, with more than 2.3 million tonnes produced worldwide in 2018, with China and Turkey as main producers [8]. China is also the largest exporter of honeys in the world, while in Europe, Portugal is the country bearing the highest number of geographical protected labels on honey [9]. Honey composition is quite variable and strictly linked to its floral source and geographical origin, but external factors, including processing, packaging and storage conditions, could play an important role. Although Italy is one of the EU countries with the highest honey production [10], the market demand for honey is higher than domestic production, therefore a substantial amount of honey is imported from elsewhere in Europe and from third-world countries, in which production does not always meet the high food safety standards required. This can lead to honey mislabeled with regard to its geographical and/or botanical origin [11].

The traceability certifying the geographical origin of food products is of primary importance for traders and producers, as well as to reinforce consumer trust. The complex task of the determination of food origin is commonly applied to control products in both customs control and self-control programs of the food industry.

Melissopalynological analysis of pollen is the most used approach for the botanical and geographical origin classification of honey, as the pollen spectrum is strictly related to the environment where the nectar is collected [12]. This analysis is often complemented by other analytical methodologies, mainly based on chromatographic techniques, to assure the honey authenticity [13]. Often, the use of conventional and targeted methods is time-consuming and not sufficient to guarantee the evaluation of complex matrices, including honeys. For this reason, the development of rapid and reliable non-targeted analytical approaches, such as fingerprinting and profiling methods, is highly demanded. Indeed, these methods combined to chemometric tools allow for the detection of a high number of metabolites, leading to samples based on their pattern.

Several analytical techniques, mainly based on nuclear magnetic resonance [5,14,15], Raman and infrared spectroscopy [16,17], mass spectrometry [18–20], electronic tongue [21,22] and electronic nose [23,24], in combination with chemometrics, have been applied to discriminate the geographical origin of honey.

The use of ambient mass spectrometry (AMS) is continuously increasing in the field of metabolomic fingerprinting as a high-throughput alternative to more traditional hyphenated methods for authentication issues [25]. Among AMS techniques, direct analysis in real-time mass spectrometry (DART-MS), being simple and requiring a very limited sample preparation, has been shown to be the most promising and versatile technique, proving to be a rapid tool in the assessment of food authenticity and food quality, also thanks to the use of fast and streamlined protocols [25–27]. Such an approach offers several advantages over the conventional techniques, including direct sample analysis in open atmosphere, high sample throughput and minimal or no sample preparation requirements, the soft ionization of a wide range of both polar and apolar compounds. Several papers have been recently published demonstrating the applicability of DART-MS to assess food authenticity and detect food adulterations of olive oil [28], beer [29], wine [30], animal fat [27], milk [31] and salmon [32]. Only one paper reported the applicability of DART-MS to the discrimination of geographical origin of food—i.e., garlic produced in Czech Republic, Spain and China [33]. Regarding honey products, DART-HRMS was used as alternative approach for the determination of 5-hydroxymethylfurfural [34,35]. To the best of our knowledge, no studies based on DART-MS have been performed to date for the assessment of geographical origin of honeys.

In this context, the aim of this study was to demonstrate the feasibility of the DART-HRMS technique for the discrimination of the geographical origin of honeys. Specifically, a rapid and suitable non-targeted DART-HRMS method in combination with multivariate statistical analysis was developed and validated to discriminate two monofloral honeys varieties (i.e., chestnut and acacia) for their geographical origin (i.e., Italy and Portugal and Italy and China, respectively). Different statistical classification models were investigated and applied to the analysis of honey samples and performance results were compared.

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

#### *2.1. Chemicals and Reagents*

Methanol (HPLC grade) was purchased from Sigma-Aldrich (Milan, Italy). Ultrapure water was produced by a Milli-Q® Direct system (Merck KGaA, Darmstadt, Germany). Helium (99.9995% purity) was provided by Sapio S.r.l. (Bari, Italy). Regenerate cellulose (RC) syringe filters with 0.2 µm of porosity were purchased by from VWR International (Milan, Italy). OpenSpot (OS) Sample Cards were purchased by Ion Sense Inc. (Saugus, MA, USA).

#### *2.2. Honey Samples*

A total of 234 monofloral honey samples commonly found in marketplaces and collected in different countries with certified origins were selected for this study. Specifically, 117 chestnut honey samples were collected from Italy (39) and Portugal (78), while 117 acacia honey samples were collected from Italy (78) and China (39). The authenticity of the monofloral honeys was assessed by internal certified protocols performed by Coop Italia Soc. Cooperativa (Casalecchio di Reno, Italy) which provided samples. Only honey samples produced in seasons 2017–2018 were taken into account.

#### *2.3. Sample Preparation*

Sampling and homogenization of honey samples were performed according to AOAC 920.180 protocol [36]. For sample preparation, a rapid protocol aimed at retaining as many honey metabolites as possible—thus to obtain most comprehensive spectra applicable for discriminating between different geographical origin—was optimized for the DART-HRMS analysis. In particular, an aliquot (1 g) of homogenized honey was added to a mixture of MeOH/H2O (1:1, *v*/*v*), (50 mL) and the sample was vortexed for 3 min. After filtration using 0.2 µm RC syringe filter, the filtered extract was directly analyzed by DART-HRMS.

#### *2.4. DART-HRMS Analysis*

DART-HRMS analyses were carried out by using a DART ionization source SI-140-GIST (DART Thermo Ion Max Vapur Interface, Ion Sense Inc., Saugus, MA, USA) coupled to an Exactive™ monostage Orbitrap™ High Resolution mass spectrometer (Thermo Fisher Scientific, San Jose, CA, USA). An aliquot (2 µL) of the honey extract was placed onto the metallic grid of the OpenSpot® sample cards and kept at 60 ◦C for 5 min to facilitate solvent evaporation before its introduction into the DART source holder. The operating conditions of the DART source were: positive ion mode; helium flow of 3.2 L/min for 1 min and heated at 250 ◦C; discharge needle voltage kept at −6 kV; grid electrode voltage set to 250 V; distance between DART exit and MS inlet set at 5 mm. The operating conditions of DART source were set by DART-SVP controller (v. 4.0.x). The main settings of the Exactive™ mass spectrometer were the following: mass scan range of 100–600 m/z; resolution set at 25,000 (FWHM at m/z 200); microscan number of 4; Automatic Gain Control (AGC) Target of 3 × 10−6 ; maximum injection time (IT) of 250 ms; capillary voltage set to 30 V; tube lens voltage set to 65 V; capillary temperature kept

at 250 ◦C. Calibrations of the MS system were periodically performed by the direct infusion ESI-MS approach of the positive ion calibrating solution, provided by the manufacturer, in order to obtain a mass accuracy lower than 5 ppm. The MS system was controlled by using the Xcalibur™ v. 2.1 software (Thermo Fisher Scientific, San Jose, CA, USA).

To carry out the subtraction of the spectral background, a blank open spot card was acquired before analyzing each sample by DART-HRMS acquiring the relevant spectrum for 30 s.

#### *2.5. Data Processing and Statistical Analysis*

In the first step of data processing, DART-HRMS spectra acquired in the time range of 30 s were averaged and then subtracted of spectral background by using the Xcalibur™ software. Successively, for each honey sample, the full list of accurate m/z ratios and peak intensities obtained was exported and processed by MetaboAnalyst 3.0 (http://www.metaboanalyst.ca/) [37,38] for peak matching and alignment with mass tolerance of 0.25, imputation of missing values (replacing missing elements by using the half of the lowest measured peak intensity) and data filtering (by using Interquartile Range approach). Successively, after pre-processing obtained by data centering, the dataset was submitted to multivariate statistical analyses performed by V-Parvus software (release 2010, http://www.parvus.unige.it, Genova, Italy).

Principal Component Analysis (PCA) was used as an unsupervised technique to evaluate the presence of outliers. Specifically, PCA was applied to each single group of monofloral honey samples of different geographical origin, observing the relevant influence plots and excluding samples identified as extreme outliers. To establish the exact number of Principal Components (PCs) to be used to build PCA models, the Non-linear Iterative Partial Least Squares (NIPALS) algorithm was applied using V-fold of 10 (cross validation process, CV = 10). PCA was also used as exploratory technique to visualize the presence of natural sample clustering between monofloral honey samples in relation to their geographical origin [39].

Afterwards, three supervised pattern recognition techniques—i.e., Linear Discriminant Analysis (LDA), Partial Least Squares Discriminant Analysis (PLS-DA) and k-nearest neighbors (k-NN) [40], were exploited to classify monofloral honey samples on the basis of their geographical origin. For this purpose, the two data matrices were randomly split in two subsets: a modelling set (containing 60 samples) and a test set (containing 57 samples). Specifically, for each monofloral honey, a modelling set, composed by 30 samples for each geographical origin, was used to build the three different statistical models. Test sets, consisting of 9 Italian and 48 Portuguese chestnut honey samples and 48 Italian and 9 Chinese acacia honey samples, were used for the validation process.

The chemometric models of PCA-LDA was built by firstly performing PCA test to reduce the number of variables that exceeded the number of objects, thus preventing model overfitting; then the selected scores were used as classification variables for LDA [41,42]. Indeed, the number of variables should not exceed (n-g)/3, where n is the number of objects and g is the number of categories. Considering that modelling sets were composed by 60 objects (number of samples) and 2 categories (number of geographical origins) the maximum number of variables should be approximately 19.

The appropriate numbers of principal components, latent variables and k values, respectively, for PCA-LDA, PLS-DA and k-NN models were established by evaluating those determining the lowest prediction error rate in cross-validation (cross-validation segments, V = 10). This parameter guarantees to improve feature variables and, at the same time, to avoid model overfitting. Model performances for PCA-LDA, PLS-DA and k-NN, expressed as percentages, were compared with reference to their recognition ability—i.e., the ability to correctly classify samples of the modelling set—prediction ability in cross-validation (CV)—i.e., the ability to correctly classify samples of a test set generated in a V-fold cross validation—and prediction ability in external validation—i.e., the ability to correctly classify samples of the test set.

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

In the present study, the real-time mass spectrometry (DART-MS) combined with chemometric analysis, was used for the first time to the discrimination of two kind of monofloral honey samples, namely chestnut and acacia, based on their geographical origin. As for chestnut, Italian and Portuguese honey samples were compared to each other, while for acacia, Italian honeys were compared with samples from China.

Figure 1 reports four representative DART-HRMS average spectra, after blank subtraction, obtained for the chestnut honey extracts of Italian (Figure 1a) and Portuguese (Figure 1b) samples and acacia honey extracts of Italian (Figure 1c) and Chinese (Figure 1c) samples. *Foods* **2020**, *9*, x FOR PEER REVIEW 5 of 12 Figure 1

**Figure 1.** *Cont.*

**Figure 1.** Representative DART-HRMS positive ion spectra acquired for the sample extracts of chestnut honeys from Italy (**a**), chestnut honeys from Portugal (**b**), acacia honeys from Italy (**c**) and acacia honeys from China (**d**). NL: Normalization level. **Figure 1.** Representative DART-HRMS positive ion spectra acquired for the sample extracts of chestnut honeys from Italy (**a**), chestnut honeys from Portugal (**b**), acacia honeys from Italy (**c**) and acacia honeys from China (**d**). NL: Normalization level.

At first, a preliminary PCA was performed on pre-processed spectra of chestnut and acacia honey samples in order to explore the presence of outlier samples. PCA score plots highlighted that, in the case of chestnut samples, seven PCs described 96.3% of total variance for samples from Italy while nine PCs described 93.0% of total variance for samples from Portugal. In the case of acacia honey samples, PCA models showed that eight PCs described 93.7% of total variance for samples from Italy while nine PCs described 91.0% of total variance for samples from China. The absence of outliers in all classes was demonstrated using influence plots where the Mahalanobis distance was plotted versus sample residual. Subsequently, an explorative PCA was performed using the entire data set to obtain an overview At first, a preliminary PCA was performed on pre-processed spectra of chestnut and acacia honey samples in order to explore the presence of outlier samples. PCA score plots highlighted that, in the case of chestnut samples, seven PCs described 96.3% of total variance for samples from Italy while nine PCs described 93.0% of total variance for samples from Portugal. In the case of acacia honey samples, PCA models showed that eight PCs described 93.7% of total variance for samples from Italy while nine PCs described 91.0% of total variance for samples from China. The absence of outliers in all classes was demonstrated using influence plots where the Mahalanobis distance was plotted versus sample residual.

of the data distribution for each monofloral honey. Figure 2 shows the PCA score plot (PC1 vs. PC2) obtained for chestnut honey samples (Figure 2a) and for acacia honey samples (Figure 2b). A discrete visual clustering of the objects on the basis of their geographical origin was observed for chestnut honeys (PC1 and PC2 explained 88.6% and 10.2% of the total variance, respectively), while overlapping regions were observed for acacia honeys with a modest clustering for their geographical origin (with 88.4% and 9.7% of the total variance explained by PC1 and PC2, respectively). Additionally, by analyzing the score plots of the remaining PCs no visual clusterization was observed. Subsequently, an explorative PCA was performed using the entire data set to obtain an overview of the data distribution for each monofloral honey. Figure 2 shows the PCA score plot (PC1 vs. PC2) obtained for chestnut honey samples (Figure 2a) and for acacia honey samples (Figure 2b). A discrete visual clustering of the objects on the basis of their geographical origin was observed for chestnut honeys (PC1 and PC2 explained 88.6% and 10.2% of the total variance, respectively), while overlapping regions were observed for acacia honeys with a modest clustering for their geographical origin (with 88.4% and 9.7% of the total variance explained by PC1 and PC2, respectively). Additionally, by analyzing the score plots of the remaining PCs no visual clusterization was observed.

**Figure 2.** *Cont.*

**Figure 2.** PC1 vs. PC2 scatter plots for monofloral chestnut (**a**) and acacia (**b**) honey samples.

These results were confirmed by analyzing the Fisher weight (FW) values of the principal components, which measure the between-class variance/within-class variance ratio. Indeed, FW values resulted to be 2.64 for the PC1 of chestnut honeys samples and lower than 1 for all the remaining PCs of both data sets (data not shown). These results indicated that the PCA was not able to discriminate honey samples on the basis of their geographical origin; therefore, it was necessary to treat data with three different supervised discriminant techniques—i.e., PCA-LDA, PLS-DA and k-NN. These classification techniques were tested on both chestnut and acacia honey samples previously split into two subsets: a modeling set and a test set. Overall, results are indicated in Tables 1 and 2, for chestnut and acacia honeys, respectively.

**Table 1.** Model performances in terms of recognition, cross validation (CV) prediction abilities and external prediction to classify chestnut honeys based on their geographical origin.

a : Italy; b : Portugal; c : Cross Validation; d : Principal Components—Linear Discriminant Analysis; e : Partial Least Squares—Discriminant Analysis; f : k-nearest neighbors.

As for LDA, PCA was used as strategy for variable reduction and to avoid model overfitting. The number of PCs (seven and nine for chestnut and acacia honeys, respectively) to be used to build the PCA-LDA models was selected on the basis of the error in prediction cross validation that has to be the lowest (CV procedure with V = 10). The PCA-LDA models provided mean value of recognition ability of 98.4% for chestnut honeys (Table 1) in both classification and CV prediction and 95.0 and 93.4% for acacia honeys (Table 2) in classification and CV prediction, respectively. The model applicability was also tested by using the test set providing mean prediction abilities of 90.3 and 89.2%, for chestnut and acacia honeys, respectively (Tables 1 and 2).

**Table 2.** Model performances in terms of recognition, cross validation (CV) prediction abilities and external prediction for all models built to classify acacia honeys based on their geographical origin.

a : Italy; b : Portugal; c : Cross Validation; d : Principal Components—Linear Discriminant Analysis; e : Partial Least Squares—Discriminant Analysis; f : k-nearest neighbors.

PLS-DA was applied as an alternative multivariate statistical approach of classification offering the advantage to avoid variables reduction processes. Specifically, by applying a 10-fold cross-validation, 10 and 12 latent variables (LVs) were found to produce the optimal model complexity for chestnut and acacia honey data sets, respectively. In these conditions, mean recognition rates were higher than 96.7% in both cases (Tables 1 and 2). Specifically, all Italian samples was correctly classified, while one Portuguese and two Chinese samples were not correctly assigned. The mean CV prediction rates were 96.7 and 95.0%, for chestnut and acacia honeys, respectively. In addition, mean prediction abilities of 89.2% and 85.8% for chestnut and acacia honeys samples, respectively, were obtained for the external validation procedure (Tables 1 and 2).

In the case of k-NN, the prediction error rate in cross-validation (V = 10) was calculated for each different k value. The smallest k value determining the lowest error was 3 for both data sets and therefore it was selected as the optimal value. The k-NN models provided mean recognition abilities in the range between 95.0–98.4%, while CV predictions were of 98.4 and 91.7%, for chestnut and acacia honeys samples, respectively. Finally, mean prediction abilities of 91.4% and 90.3% were obtained in the external validation for chestnut and acacia honeys, respectively.

The results herein obtained were in accordance with a similar study focused on the geographical authentication of Italian honey based on an NMR-metabolomic approach [5]. The authors developed a PLS2-DA model able to correctly discriminate 100% of Italian honeys from Eastern European ones. In another study, MIR analysis in combination with a PCA-LDA model were found able to distinguish geographical origins of monofloral honeys from Switzerland, Germany, and France, with prediction abilities ranged from 76 to 100% [17] although only a limited number of samples was used for the analysis.

In the current study, the DART-HRMS untargeted approach coupled with three supervised techniques, such as PCA-LDA, PLS-DA and k-NN, were investigated for discriminating Italian chestnut and acacia honey from Portuguese and Chinese samples. The results showed that all developed models provided acceptable and comparable prediction abilities, highlighting the robustness of the entire method, its applicability being unaffected by the statistical approach used to assess the authenticity of unknown samples. Moreover, these results demonstrated that DART-HRMS technique provides informative experimental data useful to build up appropriate models for the discrimination of monofloral honey samples on the basis of their geographical origin.

#### **4. Conclusions**

In this study, a rapid, easy-to-perform and low-cost method based on DART-HRMS analysis combined to multivariate statistical analysis was successfully developed and applied to classify monofloral honeys for their geographical origins, such as Italy and Portugal for chestnut samples and Italy and China for acacia samples. Specifically, three supervised approaches—i.e., PCA-LDA, PLS-DA and k-NN were evaluated. All tested models provided high and comparable recognition and prediction abilities in cross-validation and external validation, with mean values ranging from

89.2% and 98.4%. The performances of the proposed DART-HRMS method makes it an effective tool to assess the authenticity of honeys, for both industries of sector against unfair advantages of competitors and control bodies to fight food frauds. Future efforts will be directed to improve the current predictive models in order to discriminate honey samples from different production seasons and identify potential markers useful for developing a DART-HRMS target method aimed at honey authentication. Moreover, the use of the DART-HRMS approach, generating huge information in a single run, would be a useful tool for discriminating honey samples with similar organoleptic characteristics but different quality levels.

**Author Contributions:** Conceptualization: V.L., M.A., A.F.L. and L.M.; methodology, V.L. and L.M.; validation, V.L. and E.D.A.; formal analysis, G.M.F. and A.D.G.; writing—original draft preparation, V.L.; writing—review and editing, V.L., E.D.A., G.M.F., A.D.G., A.F.L., M.A. and L.M.; supervision, A.F.L. and L.M.; funding acquisition, M.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** The present research has received funding from the European Union's Seventh Framework Programme for research, technological development and demonstration, under grant agreement No. 613688 2 "Food Integrity". The equipment used in this work was funded by project BioNet – PTP – "Biodiversità per la valorizzazione e sicurezza delle produzioni alimentari tipiche pugliesi (codice n. 73, PO Regione Puglia FESR 377 2000-2006)".

**Acknowledgments:** The authors thank Salvatore Cervellieri for his support in the statistical analysis and Fernando Gottardi from Coop Italia (Casalecchio di Reno, Italy) for providing the honey samples analyzed in the present study.

**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/).

## **Machine Learning Approaches Applied to GC-FID Fatty Acid Profiles to Discriminate Wild from Farmed Salmon**

**Liliana Grazina 1 , P. J. Rodrigues 2 , Getúlio Igrejas 2 , Maria A. Nunes 1 , Isabel Mafra 1,\* , Marco Arlorio 3 , M. Beatriz P. P. Oliveira 1 and Joana S. Amaral 4,\***

Received: 23 September 2020; Accepted: 4 November 2020; Published: 7 November 2020

**Abstract:** In the last decade, there has been an increasing demand for wild-captured fish, which attains higher prices compared to farmed species, thus being prone to mislabeling practices. In this work, fatty acid composition coupled to advanced chemometrics was used to discriminate wild from farmed salmon. The lipids extracted from salmon muscles of different production methods and origins (26 wild from Canada, 25 farmed from Canada, 24 farmed from Chile and 25 farmed from Norway) were analyzed by gas chromatography with flame ionization detector (GC-FID). All the tested chemometric approaches, namely principal components analysis (PCA), t-distributed stochastic neighbor embedding (t-SNE) and seven machine learning classifiers, namely k-nearest neighbors (kNN), decision tree, support vector machine (SVM), random forest, artificial neural networks (ANN), naïve Bayes and AdaBoost, allowed for differentiation between farmed and wild salmons using the 17 features obtained from chemical analysis. PCA did not allow clear distinguishing between salmon geographical origin since farmed samples from Canada and Chile overlapped. Nevertheless, using the 17 features in the models, six out of the seven tested machine learning classifiers allowed a classification accuracy of ≥99%, with ANN, naïve Bayes, random forest, SVM and kNN presenting 100% accuracy on the test dataset. The classification models were also assayed using only the best features selected by a reduction algorithm and the best input features mapped by t-SNE. The classifier kNN provided the best discrimination results because it correctly classified all samples according to production method and origin, ultimately using only the three most important features (16:0, 18:2n6c and 20:3n3 + 20:4n6). In general, the classifiers presented good generalization with the herein proposed approach being simple and presenting the advantage of requiring only common equipment existing in most labs.

**Keywords:** authenticity; fish; *Salmo salar* L.; fatty acids; mislabeling; chemometrics; machine learning

#### **1. Introduction**

In recent decades, the consumption of fish has been increasingly recommended due to its health benefits, mainly related to the prevention of cardiovascular diseases [1]. In particular, fatty fishes from cold waters, such as salmon, are frequently rich in polyunsaturated fatty acids (PUFA), including the essential fatty acids linoleic (18:2n6) and α-linolenic (18:3n3), but also in several omega-3 PUFA such as eicosapentaenoic (EPA, C20:5n3) and docosahexaenoic (DHA, 22:6n3) acids. Besides being components of cell membranes, omega-3 PUFA are involved in the biosynthesis of eicosanoids and have been shown to influence health by affecting cell signaling cascades and gene expression, resulting in decreased expression of inflammatory and atherogenesis-related pathways [2,3]. Moreover, different studies showed that omega-3 PUFA play an important role in altering blood lipid profiles and associate their consumption with improved cardiovascular function and decreased risk of atherosclerosis and peripheral arterial disease [2,3].

In addition to these benefits, fish is largely consumed for its nutritional value and sensory aspects, making it one of the most traded food commodities. In this sense and considering that the world's wild fish stocks are limited, the production of farmed fish has been steadily increasing in recent last years. In fact, according to the Food and Agriculture Organization (FAO) Globefish Highlights, world fisheries capture was 92.5 million tonnes in 2017 with this figure expected to decrease to 91.3 million tonnes in 2019, by the contrary, fish capture arising from aquaculture is expected to grow from 80.1 to 86.5 million tonnes in the same period [4]. Concerning salmon, from 2000 to 2014, a much stronger increase was verified for aquaculture production (from 898,800 to 2,326,300 tonnes) compared to that of the world's capture of wild salmon (from 728,000 to 879,000 tonnes) [5]. Aquaculture allows wider consumer access to fish generally at more affordable costs, though it is known that fatty acid composition can significantly vary according to its production method (wild vs. aquaculture). Particularly for salmon, it has been reported that wild salmon generally present higher contents of valuable omega-3 PUFA [6–8]. This aspect, together with particular organoleptic characteristics, has driven several consumers to prefer wild salmon. Considering the limited availability of this type of salmon and its growing demand, prices have been increasing significantly, resulting in this product being prone to adulteration by origin mislabelling or even substitution with other lower-cost fish [9–11]. Whereas fish species authentication can be performed using well established and straightforward DNA-based methods [12], different approaches have been proposed so far to assess the origin of fish with respect to production method. These include, mainly, the use of nuclear magnetic resonance (NMR) [13,14], isotope ratio analysis [15,16], lipidic profile [17,18] or a combination of these [11,19–21]. Excellent discrimination (100%) between wild and farmed Atlantic salmon was reported by Aursand et al. [13] by applying support vector machines (SVM) to data obtained by 13C NMR. In another study of the same group, the lipid extract was analyzed by 13C NMR and by gas chromatography with flame ionization detector GC-FID for fatty acid composition to discriminate between wild and farmed Atlantic salmon and assign the origin of the aquaculture samples to the farms included in the study [19]. The application of chemometrics to the reference farmed fish showed very good results for both approaches, but, surprisingly, slightly better for GC-FID data. The use of stable isotope analysis based on isotope ratio mass spectrometry (IRMS) is also a promising approach, especially when combined with chemical composition analysis, notably fatty acids [11,15,21]. Yet, previous works have demonstrated that lipidic profile is sufficient to establish the production method of salmon samples, particularly when combined with chemometric analysis [8,19,20]. Recently, Fiorino et al. [8] analyzed the lipid extracts obtained from a total of 100 samples of farmed and wild salmon by direct analysis in real time (DART) coupled to high resolution mass spectrometry (HRMS). The proposed methodology showed to be fast and allowed a good discrimination between the two groups (wild vs. farmed), though without differentiating the geographical origin of the farmed fish. Moreover, the referred approach requires advanced and expensive equipment, which is not available in most control quality/analytical laboratories. In the present study, the fatty acid composition of the same samples of wild and farmed salmon used in the work of Fiorino et al. [8] was analyzed by GC-FID, an affordable equipment commonly available

in most laboratories. Subsequently, the obtained data were submitted to advanced chemometric analysis to establish the most suitable classifier able to discriminate the origin of salmon samples (wild vs. farmed, and the geographical origin among farmed samples) with the minimum possible computational effort.

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

#### *2.1. Samples*

In this study, a total of 100 authentic salmon samples obtained in the framework of the EU-funded project FOODINTEGRITY (Working Package 18) were analyzed. The samples included 26 wild salmon captured in Canada, and 74 farmed salmon samples from aquaculture farms in Canada (25), Norway (25) and Chile (24). No information was available about the gender of each specimen, neither of the diet or farming conditions used. The samples (entire fish) were transported frozen to the laboratory (Meriex Nutriscience, Chicago, IL, USA), allowed to defrost overnight at refrigerated temperature, and filleted in a cold room (4 ◦C). After removing the bones and skin, the muscles were grinded and distributed in labelled glass jars containing approximately 200 g each. The jars were immediately frozen and then shipped under freezing conditions (−20 ◦C) to the participating laboratories in different countries. After arriving, the samples were kept at −20 ◦C and submitted to lipid extraction as soon as possible.

#### *2.2. Lipid Extraction*

Lipids were extracted based on the Bligh and Dyer protocol [22] with some modifications. Briefly, about 13 g of each minced fillet were added with 13 mL of NaCl (1%) and 100 µL of butylated hydroxytoluene (BHT) (0.01% in n-hexane) to avoid oxidation, and homogenized for 1 min using an Ultra-Turrax at 13,500 rpm, keeping a low temperature by immersing the tube with the sample on ice. After that, 2.5 mL of the homogenate was transferred to a new tube and added with 2.5 mL of chloroform and 5 mL of methanol, both refrigerated. The solution was mixed vigorously by vortexing for 2 min. After centrifuging (4000 rpm, 15 min at 4 ◦C) the upper layer was discarded, and an additional 2.5 mL of refrigerated chloroform was added. After vortexing for 30 s and centrifuging under the same conditions, the chloroformic phase was transferred into a new tube and centrifuged (4000 rpm, 5 min at 4 ◦C). Finally, the chloroformic phase was collected into a previously weighted vial, flushed with a nitrogen stream and stored at −20 ◦C until further analysis. Each sample was submitted to independent extractions (*n* = 3).

#### *2.3. Fatty Acids Analysis by GC-FID*

Fatty acids were methylated using acid-catalysed trans-methylation with BF3 [23]. Firstly, the lipidic chloroformic extracts, previously stored at −20 ◦C, were dried under nitrogen and the tubes weighted to calculate the extraction yield. After dissolving the obtained lipids in 1 mL of n-hexane, for each sample, the volume containing 12.5 mg of lipids was transferred for a glass tube and dried under nitrogen. After adding 100 µL of BHT (0.01% in n-hexane) to prevent oxidation phenomena, fatty acid methyl esters were prepared. For that purpose, 1.25 mL of KOH (0.5 M) in methanol were added and the mixture was heated for 10 min at 100 ◦C after vortex-mixing vigorously. After cooling down the tubes, 1.0 mL of 14% boron-trifluoride in methanol (≥99.0% purity) (Sigma-Aldrich, Steinheim, Germany) was added to the solution, which was homogenized by vortexing, and the tubes heated again for 30 min at 100 ◦C. After completely cooling down the tubes in ice, 2.0 mL of n-hexane high performance liquid chromatography HPLC grade (Merck, Darmstadt, Germany) was added and the solution was vortex-mixed. Then, 1.0 mL of a saturated NaCl solution was added, followed by vigorous mixing and then by centrifuging for 5 min at 3000 rpm to obtain a clear upper phase. After that, 1.5 mL of supernatant was transferred to a new vial, added with anhydrous Na2SO4 and approximately 1.0 mL of FAME solution was transferred to an injection vial.

GC-FID analysis was carried out in a Shimadzu GC-2010 Plus gas chromatograph equipped with a Shimadzu AOC-20i auto-injector and a flame ionization detector (Shimadzu, Japan). FAME separation was achieved on a CP-Sil 88 silica capillary column (50 × 0.25 mm i.d., 0.20 µm, Varian, Middelburg, The Netherlands). The injector and detector temperatures were 250 ◦C and 270 ◦C, respectively. The oven parameters were set as follows: an initial temperature of 150 ◦C was increased at 3 ◦C/min to 160 ◦C and held for 2.0 min, then it was increased at 3 ◦C/min to 220 ◦C and held for 10 min. Helium was used as the carrier gas at a flow rate of 1 mL/min, and 1 µL of sample was injected using a split ratio of 1:50. Identification of compounds was performed by comparison of their retention times with those of authentic standards mixtures, namely 37 component FAME mix (certified reference material CRM47885) and PUFA nº.1 Marine source (standard 47,033) both from Supelco (Bellefonte, PA, USA). In addition, the fatty acid cis-11-octadecenoate (C18:1n7) was identified with and individual standard also purchased from Supelco. The results were expressed as the relative percentage of each fatty acid, calculated based on the chromatographic peak area. Each lipid extract was injected in duplicate.

#### *2.4. Chemometric Analysis*

#### 2.4.1. Dataset