1. Introduction
Honey is a natural sweet food produced by bees (Apis mellifera) with numerous nutraceutical and therapeutical properties [1,2,3,4,5,6,7]. Several factors, including botanical origin, environmental conditions, and bee species, impact the composition of honey and play a key role in determining its sensorial and health-promoting properties [7].
Owing to its unique qualities, honey is particularly vulnerable to unauthorized food manipulation. Adulteration and mislabeling are the most prevalent forms of honey fraud [8]. According to European (EU) legislation [9], honey labeling must indicate EU or non-EU origins, whereas botanical and local geographical origins are optional [10]. Nevertheless, because these factors strongly affect the economic value and price of honey, they are frequently falsified. This type of counterfeiting has a remarkable economic impact, especially on small beekeepers who tend to produce rare or unifloral honey [11]. For these reasons, there is considerable scientific interest in developing innovative analytical methods to verify the honey authenticity [12].
Honey’s botanical origin is generally ascertained through melissopalynological analysis. However, this method is ineffective if the honey has been filtered. On the other hand, there are no standard methods for determining geographical origin. Current methods involve advanced analytical tools coupled with chemometrics [13]. Physicochemical [14,15,16,17], elemental [18], isotopic [19], chromatographic, and hyphenated mass spectrometry methods [20,21], NMR [22,23], DNA-based [24,25], electronic sensing [26,27], and spectroscopic [28,29,30,31,32] techniques are nowadays the most used for predicting botanical and geographical origin, or detecting adulterants [33]. Usually, botanical issues are investigated using chromatographic techniques and physicochemical analyses [12,14,34]. However, while vibrational spectroscopy is commonly used to detect adulterated honeys [12,33], it is also highly effective in predicting the authenticity of almost all types of honey. Furthermore, the use of portable instrumentation enhances the versatility of this approach [35,36]. Finally, elemental and isotopic methods are prevalently used for geographical origin discrimination [19,37], but they are also useful in determining botanical origin [37,38]. Elements in honey reflect the soil composition [39], flora [40], and anthropogenic activities [41]. The translocation of elements is influenced by season [42], bee species [43], and, in general, by environmental conditions [44]. Furthermore, honey reflects the urbanization grade [39,45,46]. All these factors determine the elemental fingerprint of honey, allowing its authentication in terms of both geographical and botanical origins [47].
The elemental composition of honey has been characterized worldwide. Recent studies investigated honeys from Argentina [48], Bulgaria [40,49,50], China [43,51,52,53], Croatia [42,49,54,55], Denmark [56], Egypt [57], Eritrea [58], Ethiopia [59], France [60], Germany [56], Greece [37,57,61], Hungary [62,63], Israel [64], Italy [37,38,46,65,66,67,68,69,70,71,72,73], Kazakhstan [73], Kosovo [74], Montenegro [75], Morocco [49,57], New Zealand [45,56], Poland [37], Romania [37,60], Serbia [74], Slovenia [49], Spain [57,76,77,78,79], Thailand [39], Tunisia [68], and Turkey [16,49]. Typically, data management and elaboration involve multivariate analysis and machine-learning techniques, including Principal Component Analysis (PCA) [38,59], Cluster Analysis (CA) [37,39], Linear Discriminant Analysis (LDA) [57,59], Partial Least Squares regression (PLS) [37,52], Decision Tree Analysis (DTA) [56], Self-Organizing Maps (SOMs) [40], and Soft Independent Modelling of Class Analogy (SIMCA) [60]. Few papers in the literature report cases where elemental fingerprinting has discriminated honeys according to their botanical and geographical origin [37,42,49,60]. Research often focuses on authenticating honeys from various origins based on botanical sources or regions with unique climates, soils, or levels of urbanization. European legislation requires differentiation between EU and non-EU honeys, which is usually feasible due to the numerous environmental factors that affect their elemental fingerprints. On the other hand, discriminating honey from nearby areas is a more challenging task.
The main goal of this study is to distinguish unifloral and multifloral honeys from two nearby regions. In this context, Sardinia (Italy) and Spain have been chosen as case studies. Geologically, the Corsica–Sardinian microplate separated from the Iberian Peninsula during the Miocene. Therefore, the soil composition of Sardinia is more similar to that of southeastern Spain than Italy [80]. In addition, the two regions have similar climatic conditions.
A comprehensive selection of multifloral and unifloral honeys from common botanical species, such as rosemary and eucalyptus, was considered for this purpose. The physicochemical characterizations of eucalyptus and rosemary honeys from Spain [81,82] and Italy [83,84] have been reported in literature. Possible differences among the eucalyptus and rosemary honeys are related to their diastase activity and acidity parameters. The diastase activity of Italian eucalyptus honey samples is higher than that of Spanish ones, whereas the acidity parameters of Spanish eucalyptus honey samples are higher than those of Italian ones. Additionally, Spanish rosemary honey samples exhibit higher diastase activity than Italian ones. Furthermore, some typical unifloral honeys from Sardinia, such as strawberry tree, asphodel, and thistle honey [85], were analyzed. Strawberry-tree honey is famous worldwide for its unique “bitter” taste [86] and for its healing properties [87]. The chemical composition of strawberry-tree honey confirms its botanical origin due to the presence of high concentrations (several hundred mg kg−1) of homogentisic acid [86,88], a compound not found in other unifloral honeys. The botanical origin of asphodel honey is also guaranteed by the peculiar presence of methyl syringate in it [89]. Among all the Sardinian unifloral honeys considered in this study, the thistle honey is characterized by the lowest values of pH, electrical conductivity, and color [14]. The honeys’ elemental fingerprints were determined using a validated Inductively Coupled Plasma Mass Spectrometry (ICP-MS) method [38]. Four major elements (Na, Mg, K, Ca), twenty-three among trace and toxic elements (Ag, As, Ba, Be, Bi, Cd, Co, Cr, Cu, Fe, Hg, Li, Mn, Mo, Ni, Pb, Sb, Sn, Sr, Te, Tl, V, Zn) and fourteen lanthanides (La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, Lu) were analyzed. The data were processed using PCA for multivariate data visualization, whereas LDA and Random Forest were used for honey classification.
2. Materials and Methods
2.1. Honey Samples
The study analyzed honeys from two areas in the Western Mediterranean: Spain (SPA) and Sardinia (Italy, ITA). The geographical and botanical origins of the honey are shown in Figure 1.
In total, 247 honey samples from Spain (SPA = 73) and Sardinia (ITA = 174) were examined. The sample set consisted of both multifloral and unifloral honeys. Among the common honeys from the two regions, multifloral (MUL, SPA = 34, ITA = 35), eucalyptus (EUC, SPA = 13, ITA = 30), and rosemary (ROS, SPA = 26, ITA = 6) were included. Additionally, three characteristic unifloral Sardinian honeys, asphodel (ASP, ITA = 33), strawberry tree (STR, ITA = 31), and thistle (THI, ITA = 39), were considered for the geographical attribution. The botanical origin of the samples was determined by melissopalynological analysis. The collection was recorded between 2020 and 2022, reflecting the flowering and seasonality of the botanical sources. In general, eucalyptus, thistle, and multifloral honeys were produced in spring, summer, and fall. Rosemary and strawberry-tree honeys were collected in fall and winter, while asphodel honeys were produced in winter and spring. Sardinian honeys were gathered throughout the island (Figure 1), whereas the Spanish honeys were from Andalusia, Aragon, Asturias, Cantabria, Castilla la Mancha, Castilla-Leon, Catalonia, Extremadura, Balearic Islands, Navarre, and Basque Country (Figure 1). Previous studies [81] have reported details on the chemical–physical and sensory characteristics of the honeys, such as color, pH, moisture, taste, and botanical markers. Honeys were stored in the dark at 4 °C until analysis.
2.2. Instrumentation and Reagents
The elemental analysis was performed using a NexION 300X ICP-MS spectrometer from Perkin Elmer (Milan, Italy). The spectrometer was equipped with an S10 autosampler, glass concentric nebulizer, glass cyclonic spray chamber, and kinetic energy discrimination (KED) collision cell. Microwave acid digestion was performed using an ultraWAVE™ from Milestone (Sorisole, Italy) equipped with a single reaction chamber (SCR) system, a fifteen-position rotor, and polytetrafluoroethylene (PTFE) vessels (15 cm3). Dry ashing was performed using a Controller P320 muffle from Nabertherm (Lilienthal, Germany). Samples were homogenized before analysis using an Ultraturrax mixer model T18 (Staufen, Germany). Nylon filters (pore diameter, 0.22 μm), syringes, metal-free polypropylene tubes (15 and 50 cm3), and porcelain crucibles (150 cm3) were supplied by VWR (Milan, Italy).
A MilliQ Plus System (Millipore, Vimodrone, Italy) was used to produce type I water (resistivity > 18 MΩ cm−1). Nitric acid (67–69% w/w, NORMATON® for ultra-trace analysis) and hydrogen peroxide (30% w/w, NORMATON® for ultra-trace metal analysis) were supplied by VWR (Milan, Italy). Periodic table mix 1 for ICP (TraceCert®, 33 elements, 10 mg dm−3 in 10% HNO3), periodic table mix 3 for ICP (TraceCert®, 16 elements, 10 mg dm−3 in 5% HNO3), Rh solution (1000 mg dm−3 in 3% HNO3), apple leaves NIST SRM® 1515 and mussel tissue BCR 668 were from Sigma-Aldrich (St. Louis, MO, USA). Single standard solutions of Na, Mg, K, Ca, Hg, Mo, Sb, and Sn (100–1000 mg dm−3 in 2–5% HNO3) were obtained from Carlo Erba (Milan, Italy).
2.3. Sample Preparation
Sample preparation involved the use of microwave acid digestion and dry ashing techniques. Microwave acid digestion was utilized to determine macroelements, trace elements, and toxic elements, while dry ashing was employed for lanthanide analysis [70]. Microwave acid digestion, performed according to a previously described method [38], was optimized by a 22 full factorial experimental design. In this manner, the residual carbon and acidity levels were minimized, preventing the need for unnecessary dilutions and reducing the matrix effect. Initially, the samples were homogenized at 40 °C. Then, approximately 0.700 g of honey was weighed in 15 cm3 PTFE vessels and treated with 0.5 cm3 of HNO3, 3 cm3 of H2O2, and 4 cm3 of type I water. After digestion at 240 °C, samples were collected, diluted to 15 cm3, and filtered before analysis. Dry ashing was performed weighing about 5.0 g of honey in porcelain crucibles (150 cm3). After ashing at 600 °C, samples were treated with 10 cm3 of 5% HNO3 aqueous solution, diluted to 15 cm3, and filtered before analysis. Table S1 reports the operational conditions of both methods.
2.4. Elemental Analysis
A previously developed and validated procedure [38] was used on a NexION 300X ICP-MS (Perkin Elmer) to perform elemental analysis. Here, detailed information regarding instrumental parameters, elemental settings, method assessment, performance, quality control, and validation is reported. The literature method for the analysis of trace and toxic elements (i.e., Ag, As, Ba, Be, Bi, Cd, Co, Cr, Cu, Fe, Hg, Li, Mn, Mo, Ni, Pb, Sb, Sn, Sr, Te, Tl, V, Zn) in honey [38] was implemented to analyze macro elements (i.e., Na, Mg, K, Ca) and lanthanides (i.e., La, Ce, Pr, Nd, Sm, Eu, Gd, Dy, Ho, Er, Tm, Yb, Lu). Tables S2 and S3 report the instrumental conditions and elemental settings for the analysis of macro elements and lanthanides, respectively. Trueness was evaluated by analyzing two certified reference materials, apple leaves NIST SRM® 1515, and mussel tissue BCR 668. The results are presented in Table S4.
2.5. Statistical Analysis
Data analysis was conducted using the statistical freeware software R (v. 4.3.1) run in the free integrated development environment R-Studio (v. 2023.3.1), GraphPad Prism (v. 9.1.0 221), and Chemometric Agile Tool (CA) [90]. For data visualization, PCA was performed removing all elements that were rarely quantified. Compositional data analysis (CoDa) was applied as a data pre-treatment before PCA to improve the interpretability using centered log-ratio transformation [91]. To remove missing values (i.e., those below the limit of detection or quantification), the dl23 method (two-thirds of the limit of detection) was used [91]. To classify honeys according to botanical origin, geographical origin, and their combination, LDA was performed. A RF machine-learning algorithm was also used (R package randomForest, [92]). Data for each factor (i.e., geographical, botanical origin, and their combination) were extracted and divided into train and test sets. The train set included an equal number of samples from each group, equivalent to half of the least populated group. The remaining data constituted the test set. The classification algorithms were applied by removing chemical elements with missing (i.e., below detection or quantification limits) or repeated values, iterating 100 times both train and test sampling to increase the statistical significance. Accuracy was measured by averaging over the accuracies obtained in the 100 replicas. In particular, RF analysis was performed using 1000 trees and setting the mtry parameter equal to the square root of the number of predictors [92]. To find the smallest set of chemical elements required to efficiently discriminate between groups, the importance for the classification of variables was first measured using the mean decrease in the Gini index [92,93]. Predictors with the lower Gini index were removed and the analysis was iterated with the remaining predictors. The smallest set of chemical elements was identified as the one preceding a visible reduction in the overall classification accuracy. Statistical significance was set at p
4. Discussion
Elemental analysis allows the content of toxic and nutritional elements to be evaluated. The levels of harmful elements are similar to or frequently lower than those found in other Spanish [78,79] and Italian [46,69,72] honeys analyzed in previous studies. On the other hand, honeys have a relatively high content of minerals such as Na, Mg, K, Ca, Zn, Mn, and Cu. However, assuming a daily honey intake of 20 g, elements of nutritional interest do not cover daily requirements, while the toxic elements do not pose any health risk.
Elemental fingerprints were tested for classifying honeys according to geographical and botanical origins. PCA was performed using two different data sets pre-treatment, centered log-ratio transformation (Figure 2), and standardization (Figure S1). As previously reported [91], the centered log-ratio transformation improves the data interpretability by distributing a part of the variance explained by the first component into the other components. Consequently, loadings and scores in the first two components are more dispersed (Figure 2), allowing a more comprehensive differentiation of honey categories. The PCA results show that botanical information predominates over geographical information (compare Figure 2B with Figure 2C). Except for rosemary honeys (for which Sardinia is underrepresented), multifloral honeys tend to overlap in the graph (Figure 2E), while eucalyptus honeys are more distinguished (Figure 2D). It is hypothesized that this difference is attributable to the origin of Spanish eucalyptus honeys. These are produced mainly in the northern regions, which are geographically more different from Sardinia. Furthermore, the results of the PCA analysis indicate that geographical discrimination is facilitated, as expected, when considering different botanical origins (compare Figure 2D with Figure 2H). Finally, the score plots also allow the comparison of the honeys’ elemental compositions. Broadly, the percentage of lanthanides is lower in Spanish samples compared with Sardinian unifloral varieties. On the other hand, macroelements and trace elements are relatively less abundant in unifloral honey than in multifloral ones (Figure 2).
Regarding honey authentication, the RF algorithm performs better than LDA in classifying honeys according to their origins. Overall, the classification by geographical origin is the most accurate (90%, predictors: Na, Mg, Mn, Sr, Zn, Ce, Nd, Eu, Tb). The accuracy tends to decrease when classifying honey by botanical origin (73%, predictors: Na, K, Ca, Mn, Sr, Ce, Eu, Lu) and when combining geographical and botanical origins (65%, predictors: Na, K, Mn, Sr, Ce, Eu, Lu). Based on these results, it is possible that geographical origin may have a greater influence than botanical origin on the elemental fingerprint. However, predicting the botanical origin is challenging, and combining the two factors reduces the classification model’s accuracy.
To our knowledge, this is the first case where RF combined with elemental fingerprinting has been used to distinguish the origin of honey. Thus, the accuracy of the models cannot be easily assessed if compared with data presented in the literature. Few studies have investigated the authentication of honey in terms of both geographical and botanical origin [37,42,49,60]. Often the data were not processed with classification algorithms [42,49]. Drivelos et al. achieved excellent results in geographical classification [37]. Magdas et al. also achieved excellent results, but their classification models were obtained using isotopic markers in addition to elements [60].
In general, as expected, the models obtained indicate that unifloral varieties are easier to accurately classify than multifloral varieties of honey. Regarding elemental predictors, Na, Mn, Sr, Ce, and Eu are common to all origin classifications, while Mg, Zn, Nd, and Tb are useful for geographical classification. On the other hand, K, Ca, and Lu are relevant for botanical origin. These findings can be compared and align with those of Pavlin et al. [49], who analyzed different varieties of unifloral and multifloral honeys from Slovenia, Croatia, Bulgaria, Turkey, and Morocco [49]. They reported Mn, K, and Ca as botanical markers, and Na, Mg, and Fe as geographical ones. However, the labeling of elements as markers for the determination of a specific origin may vary depending on the honey’s origin or variety.
Regarding lanthanides, they have primarily been reported as indicators of geographical origin [19,60]. In this study, Ce, Nd, Eu, and Tb are significant for geographical classification, while Ce, Eu, and Lu are significant for botanical classification. Previously, Squadrone et al. [66] reported that lanthanides can help in discriminating unifloral and multifloral honey using the ratio of light and heavy rare earths (LREE/HREE), while Gulino et al. [70] reported that the fractionation of heavy lanthanides partially helps in geographical classification. The results of this research suggest that Ce and Eu are among the most important lanthanides in all models. Eu exhibits anomalous behavior in Spanish honeys, which is not observed in those from Sardinia. Gulino et al. [70] suggested that these anomalies could be attributed to Ba interferences during the analyses. Although it may be possible, any potential bias would affect all samples and should be highly correlated with Ba. However, the level of correlation between these parameters is low (r = 0.16), so this hypothesis can be ruled out. Based on the results obtained, Eu may be considered a reliable marker in all classification models.
