Multi-Elemental Analysis as a Tool to Ascertain the Safety and the Origin of Beehive Products: Development, Validation, and Application of an ICP-MS Method on Four Unifloral Honeys Produced in Sardinia, Italy

Despite unifloral honeys from Sardinia, Italy, being appreciated worldwide for their peculiar organoleptic features, their elemental signature has only partly been investigated. Hence, the principal aim of this study was to measure the concentration 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, and Zn) in four unifloral honeys produced in Sardinia. For this purpose, an original ICP-MS method was developed, fully validated, and applied on unifloral honeys from asphodel, eucalyptus, strawberry tree, and thistle. Particular attention was paid to the method’s development: factorial design was applied for the optimization of the acid microwave digestion, whereas the instrumental parameters were tuned to minimize the polyatomic interferences. Most of the analytes’ concentration ranged between the relevant LoDs and few mg kg−1, while toxic elements were present in negligible amounts. The elemental signatures of asphodel and thistle honeys were measured for the first time, whereas those of eucalyptus and strawberry tree honeys suggested a geographical differentiation if compared with the literature. Chemometric analysis allowed for the botanical discrimination of honeys through their elemental signature, whereas linear discriminant analysis provided an accuracy level of 87.1%.


Introduction
Honey is an ancient natural and functional food used for centuries in the traditional medicine of many cultures. Beyond monosaccharides, such as glucose and fructose, which it is composed of up to 80% w/w, honey is a very complex matrix containing many bioactive compounds such as proteins, amino acids and enzymes, organic acids, polyphenols, flavonoids, vitamins, and inorganic elements [1]. This combination makes honey an outstanding functional food with many health-promoting activities and antimicrobial, antiviral, antifungal, anticancer, and antidiabetic properties [2][3][4][5][6].
Sardinia, Italy, is the second-largest island in the Mediterranean Sea. Its extension and distance from the European and African continental shelves characterize the endemic flora and fauna. Furthermore, due to the paucity of population and industries, the poor anthropogenic pressure is considered by consumers as an intrinsic guarantee of the quality as well as origin [56,57]. Hence, the principal aim of this research was to ascertain the elemental signature of the most renowned unifloral honeys of Sardinia as a reliable tool to ensure their healthiness and guarantee their origin. Therefore, an original ICP-MS method able to simultaneously measure the total amount of 23 elements of potential health concern (i.e., Ag, As, Ba, Be, Bi, Cd, Co, Cr, Cu, Fe, Hg, Li, Mn, Mo, Ni, Pb, Sb, Sn, Sr, Te, Tl, V, and Zn) was developed and thoroughly validated. Finally, the proposed method was applied to a large sampling of honeys from asphodel, eucalyptus, strawberry tree, and thistle.

Sample Pre-Treatment
To avoid the matrix effect and maximize the signal-to-noise ratio, great attention was paid to the optimization of the sample pre-treatment step [44][45][46][47]. Although the sample amounts and the instrumental conditions depend on the specific features of the microwave system, according to the recent method proposed by Astolfi et al. [42], HNO 3 and H 2 O 2 were chosen as the acidic/oxidizing mixture due the fact of their proven capability to minimize the formation of the interfering polyatomic ions in the plasma and their efficiency in organic matter decomposition. According to Muller et al. [45], a simply 2 2 full factorial design was applied to improve the digestion efficiency and the composition of the oxidizing mixture. The optimization allowed for a reduction in the HNO 3 volume required for sample digestion, increasing the amount of the greenest hydrogen peroxide and the mass of sample digested, which is an important achievement for trace element analysis. These results were also achieved thanks to the versatile performance of the ultraWAVE digestion system. The optimized method is reported in Table 1. First, honey samples were heated until 40 • C and then homogenized by an Ultraturrax. Next, approximately 0.7 g of honey, exactly weighted on an analytical balance (±0.0001 g uncertainty), was treated with 0.5 cm 3 of HNO 3 , 3 cm 3 of H 2 O 2 , and 4 cm 3 of type I water inside a PTFE vessel of the SCR mineralization system. After the digestion, samples were diluted to a final volume of 15 cm 3 and filtered through a 0.22 µm nylon filter. Finally, they were stored at 4 • C in the dark until the analysis. The typical amount of residual acidity found in digested samples was 0.2 mol dm −3 , whereas the efficiency of organic matter decomposition (EOMD%) was generally higher than 94%.

Validation of the ICP-MS Method
The validation of the ICP-MS method proposed was accomplished in terms of the limit of detection (LoD), limit of quantification (LoQ), linearity, precision, and trueness. Table 2 reports the validation parameters of the ICP-MS method for the determination of the total amount of 23 trace elements in honey. The LoDs and LoQs of the method were calculated according to Currie [58] on 30 measurements of method blanks obtained in different analytical sessions. The LoDs were generally below 10 µg dm −3 for all of the elements except for Ba, Cu, Fe, and Zn, which were between 10 and 40 µg dm −3 .
Because of the large variability of the analyte concentrations, great attention was paid to the instrument calibration phase. Although the linearity of the ICP-MS could be extended over several orders of magnitude of concentration, in this case, this parameter was explored for each element only within the relevant operative range of concentration. For each analyte, the calibration function was the result of a linear regression of six standard solutions, three for each extreme point of the calibration range. To evaluate the experimental error and verify the linearity of the calibration function, three different solutions at an analyte concentration equal to the central point of the calibration range were analyzed. For each analyte, the difference between the experimental and predicted values was not significant (t-test, α = 0.05). Hence, keeping in consideration the random distribution of the residuals around the mean value as well as the very high coefficient of determination (R 2 always above 0.999), the linearity of the calibration function was successfully ascertained.
Precision, measured as the coefficient of variation (CV%), was assessed in terms of repeatability and intermediate precision. For this purpose, a batch of honey samples was replicated in the same way and in different analytical sessions. Table 2 reports the CV%s calculated at the analyte's average concentration. For elements in which the concentration in the honey was below the LoD, precision was measured by spiking the samples with standard solutions of analyte. Repeatability exhibited CV% between 1% (Ba) and 12% (Zn), while intermediate precision ranged between 3% (Ba and Mn) and 21% (Te). All precision parameters were in the range of CV% defined by the Horwitz's theory [59]; hence, the overall precision level of the method was acceptable.
Due to the lack of certified reference material (CRM) for the determination of trace elements in honey, trueness was evaluated by recovery tests. Hence, three aliquots of each sample were spiked with increasing amounts of a standard solution containing all analytes. All aliquots then underwent the whole analytical procedure described in Section 2.1. Furthermore, an additional aliquot for each honey sample was analyzed without any spiking. Recoveries ranged between 85% (Bi) and 130% (Hg). Quantitative recoveries (t-test, α = 0.05) were observed for most of the elements, except for As, Li, Mo, Pb, and V (underestimation bias) and for Be, Cd, and Sb (overestimation bias). Nevertheless, the recoveries obtained in these cases were also acceptable according to the AOAC guidelines [60]. In summary, the proposed method was fully and successfully validated, and the parameters here considered were comparable or better than those reported in previous works [31,35,36,41,42]. Table 3 summarizes, in terms of both the mean value and range, the amount of toxic and trace elements belonging to the four unifloral Sardinian honeys, while the full data set is available in the Supplementary Materials (Table S1). Each sample was analyzed twice. Italics represent data below the LoD; underlined values represent data below the LoQ. Table 4 shows the results of the amounts of toxic and trace elements reported for unifloral eucalyptus honeys produced in different countries. The determination of As, Cd, Cr, Cu, Fe, Mn, Ni, Pb, and Zn was performed in all studies, and the amounts of these elements in honeys from different geographical origins normally spanned over one or two orders of magnitude. The different concentrations of the most abundant elements (i.e., Cr, Cu, Fe, Ni, Pb, and Zn, significantly higher in the samples from Tunisia, or the very high concentration of Mn in the Argentinian honey) allowed to envisage the feasibility of a discrimination of these honeys according to geographical origin. However, due to the paucity of the samples analyzed in the literature [36][37][38][39], no definitive conclusion could be obtained by this comparison.

Honey Analysis
The comparison of the data set obtained from the determination of 23 elements in nine samples of strawberry tree from Croatia [40] was more reliable. In this case, 17 out of 23 elements were measured in both countries. Meaningful differences in the average concentrations of many trace elements (i.e., Co, Cu, Mo, Ni, Pb, Sb, Sn, and V) appeared to support the possibility to achieve a differentiation according to the different geographical origin. The average distributions reported in Figure 1 highlights the differences between Croatian and Sardinian strawberry tree honeys. 23 elements were measured in both countries. Meaningful differences in the average concentrations of many trace elements (i.e., Co, Cu, Mo, Ni, Pb, Sb, Sn, and V) appeared to support the possibility to achieve a differentiation according to the different geographical origin. The average distributions reported in Figure 1 highlights the differences between Croatian and Sardinian strawberry tree honeys.

Chemometric Analysis
In order to evaluate the suitability of trace elements to ascertain the origin of Sardinian honeys, principal component analysis (PCA) and linear discriminant analysis (LDA) were performed for exploratory and classification purposes, respectively.
The original data set, constituted by 133 samples and 18 variables (Ag, Be, Hg, Sb, and Te were removed because they were almost never quantified) were randomly divided into training (85 samples) and validation (48 samples) sets for internal validation.
The logarithmic transformation (log10) was used to reduce the skewness of the probability density distribution present in the original data. However, the information embodied in the loadings has changed and so it could lead to dangerous misunderstandings [61].
From the loading plot (Figure 2a), PC1 explained mainly the global concentration of the trace elements characterized by positive loadings on PC1 (Fe, V, Mo, and Mn). On the other hand, PC2 explained the contrast between the two clusters of correlated trace elements: Li, Sr, and Ba at positive loadings and Co, Cd, Cu, Zn, and Ni at negative loadings. Couples of elements, such as Mn and V, Cd and Co, and, mainly, Cu and Zn, were strongly correlated among them. Finally, Sn, Cr, Pb, and Bi were the less significant variables, as suggested by their overall low expression in both PC1 and PC2. Looking at the score plot in the plane PC1-PC2 (Figure 2b), with the objects colored according to the different botanical origin, it was noticed that the four classes were generally characterized by a specific location on the plane: asphodel honey samples (1) can be found at negative scores of both PC1 and PC2; eucalyptus samples (2) at positive scores of both PC1 and PC2; strawberry tree samples (3) at negative score of PC1 and positive scores of PC2 and, finally, thistle samples (4) can be found at positive scores of PC1 and negative scores of PC2. The information on the within-class variability of the four classes was smoothed by the logarithmic transformation, as well as the correlations between the variables. The high variability can be confirmed from the width of the ranges reported in Table 3. Furthermore, an analysis of the correlations was performed before and after the logarithmic transformation (Table S2) to ensure a correct interpretation of the variables. The coefficients highlight a strong correlation between Fe-Mn-V and Cu-Ni-Zn, confirmed by the loading plot, in

Chemometric Analysis
In order to evaluate the suitability of trace elements to ascertain the origin of Sardinian honeys, principal component analysis (PCA) and linear discriminant analysis (LDA) were performed for exploratory and classification purposes, respectively.
The original data set, constituted by 133 samples and 18 variables (Ag, Be, Hg, Sb, and Te were removed because they were almost never quantified) were randomly divided into training (85 samples) and validation (48 samples) sets for internal validation.
The logarithmic transformation (log10) was used to reduce the skewness of the probability density distribution present in the original data. However, the information embodied in the loadings has changed and so it could lead to dangerous misunderstandings [61].
From the loading plot (Figure 2a), PC1 explained mainly the global concentration of the trace elements characterized by positive loadings on PC1 (Fe, V, Mo, and Mn). On the other hand, PC2 explained the contrast between the two clusters of correlated trace elements: Li, Sr, and Ba at positive loadings and Co, Cd, Cu, Zn, and Ni at negative loadings. Couples of elements, such as Mn and V, Cd and Co, and, mainly, Cu and Zn, were strongly correlated among them. Finally, Sn, Cr, Pb, and Bi were the less significant variables, as suggested by their overall low expression in both PC1 and PC2. Looking at the score plot in the plane PC1-PC2 (Figure 2b), with the objects colored according to the different botanical origin, it was noticed that the four classes were generally characterized by a specific location on the plane: asphodel honey samples (1) can be found at negative scores of both PC1 and PC2; eucalyptus samples (2) at positive scores of both PC1 and PC2; strawberry tree samples (3) at negative score of PC1 and positive scores of PC2 and, finally, thistle samples (4) can be found at positive scores of PC1 and negative scores of PC2. The information on the within-class variability of the four classes was smoothed by the logarithmic transformation, as well as the correlations between the variables. The high variability can be confirmed from the width of the ranges reported in Table 3. Furthermore, an analysis of the correlations was performed before and after the logarithmic transformation (Table S2) to ensure a correct interpretation of the variables. The coefficients highlight a strong correlation between Fe-Mn-V and Cu-Ni-Zn, confirmed by the loading plot, in contrast, the correlations between Ba-Sr-Li and Cd-Co were emphasized by the data pre-treatment. contrast, the correlations between Ba-Sr-Li and Cd-Co were emphasized by the data pretreatment. For classification purposes, linear discriminant analysis (LDA) was used. This method is based on the description of data by means of probability density distributions under two hypotheses: (i) probability distributions are multivariate normal within all the classes; (ii) dispersion and correlation structures are the same within all the classes. Fortunately, the method is quite robust against slight deviations from these hypotheses [62].
In this case, after the pre-processing, it was possible to use LDA because the classes were similar in size and orientation and, therefore, it was possible to assume that they had similar variance and covariance matrix.
The results obtained in cross-validation and prediction are reported in Table 5. The botanical origin with the lowest percentage of correct prediction was asphodel (61.9%), whose samples were mainly confused with thistle (n = 3) and strawberry tree (n = 5), while the best results were obtained for eucalyptus (84.2%). A permutation test was performed to compare the performance obtained with the correct assignment with a series of random permutation of the classes (in this case, 100 permutations were calculated). Figure S3 shows how the model works, since in no event in the classification were the results obtained randomly. In prediction, a correct classification percentage of 87.1% was achieved, quite in agreement with the results obtained in the CV. By observing the confusion matrix (Table 6), it is possible to look at the samples characterized by the wrong assignment, highlighting that asphodel and thistle were relatively less accurate but still acceptable: strawberry tree was twice erroneously classified as thistle and once as asphodel; eucalyptus was erroneously re-classified as both thistle and asphodel; a sample of asphodel was wrongly classified as a thistle.  For classification purposes, linear discriminant analysis (LDA) was used. This method is based on the description of data by means of probability density distributions under two hypotheses: (i) probability distributions are multivariate normal within all the classes; (ii) dispersion and correlation structures are the same within all the classes. Fortunately, the method is quite robust against slight deviations from these hypotheses [62].
In this case, after the pre-processing, it was possible to use LDA because the classes were similar in size and orientation and, therefore, it was possible to assume that they had similar variance and covariance matrix.
The results obtained in cross-validation and prediction are reported in Table 5. The botanical origin with the lowest percentage of correct prediction was asphodel (61.9%), whose samples were mainly confused with thistle (n = 3) and strawberry tree (n = 5), while the best results were obtained for eucalyptus (84.2%). A permutation test was performed to compare the performance obtained with the correct assignment with a series of random permutation of the classes (in this case, 100 permutations were calculated). Figure S3 shows how the model works, since in no event in the classification were the results obtained randomly. In prediction, a correct classification percentage of 87.1% was achieved, quite in agreement with the results obtained in the CV. By observing the confusion matrix (Table 6), it is possible to look at the samples characterized by the wrong assignment, highlighting that asphodel and thistle were relatively less accurate but still acceptable: strawberry tree was twice erroneously classified as thistle and once as asphodel; eucalyptus was erroneously re-classified as both thistle and asphodel; a sample of asphodel was wrongly classified as a thistle.  The level of discrimination achieved in this contribution is less than those obtained by this research group using chemometric approaches based on untargeted physical and chemical data [49] or FT-ATR spectroscopic data [50].
Several factors could affect the ascertainment of the botanical origin. First, Sardinian unifloral honeys are typically underrepresented from a melissopalynological point of view since the pollen spectrum is often largely overlapping due to both accompanying and rare pollens [8]. This highlights a partly common botanical origin, especially for spring productions of asphodel and thistle and, secondly, also for summer and autumn production of eucalyptus and strawberry tree, respectively. Finally, a last possible factor, could be related to the geographical origin of some honeys, especially for those coming from the mining areas (i.e., central and southwestern Sardinia), where the natural characteristics of soil composition may have influenced the concentration of some elements.
In addition, a second LDA was unsuccessfully used in the attempt to achieve a geographical classification of Sardinian samples belonging to the same floral origin (data not reported). Probably, the reason for this failure was related both to the difficulty to localize the exact position of the hives [8] and to the extreme geochemical and pedological variability of the soil of Sardinia [63].
Summarizing, the elemental signature provided rather robust descriptors for the botanical discrimination of Sardinian unifloral honeys, whereas the intra-regional geographical discrimination was currently compromised and should be investigated by acquiring additional data. Hence, efforts will be paid to carefully solve these issues. Moreover, additional samples of both eucalyptus and strawberry tree honeys produced outside Sardinia will be gathered in order to achieve a reliable geographical discrimination.

Honey Samples
The honeys were mainly gathered in the last year directly from local beekeepers and, in a minor amount, purchased in Sardinian markets. The collection, summarized in Figure 3, consisted of 133 samples of the four most renowned unifloral varieties of honey: asphodel (n = 33), eucalyptus (n = 30), strawberry tree (n = 31), and thistle (n = 39). Samples were stored in the dark at 4 • C until the analysis. The botanical origin of each sample, primarily based on information directly provided by beekeepers, has been confirmed by the melissopalynological analysis, which provided data within the relevant ranges measured by Floris et al. [8].
Briefly, Sardinia is mainly hilly, therefore, 60% (n = 80) of the honeys are from the hills, 36% (n = 48) from the plains, and only 4% (n = 5) from the mountains. The production areas are divided into rural areas 50% (low urbanization degree, n = 67), slightly urbanized areas 45% (medium urbanization degree, n = 62) and urbanized 3% (high urbanization degree, n = 4) [64]. The odd distribution of the samples in the region, as reported in Figure 3, is representative of the different production areas and the distribution of the botanical sources. Briefly, Sardinia is mainly hilly, therefore, 60% (n = 80) of the honeys are from the hills, 36% (n = 48) from the plains, and only 4% (n = 5) from the mountains. The production areas are divided into rural areas 50% (low urbanization degree, n = 67), slightly urbanized areas 45% (medium urbanization degree, n = 62) and urbanized 3% (high urbanization degree, n = 4) [64]. The odd distribution of the samples in the region, as reported in Figure  3, is representative of the different production areas and the distribution of the botanical sources.

Reagents and Standard Solutions
In all the analytical phases, type I water (resistivity > 18 MΩ cm −1 ), produced by means of a MilliQ plus System (Millipore, Vimodrone, Italy) was used. Nitric acid (67-69% w/w, NORMATON ® for ultra-trace analysis) and hydrogen peroxide (30% w/w, NOR-MATON ® for ultra-trace metal analysis) were from VWR (Milan, Italy). The multi-element standard periodic table mix 1 for ICP (TraceCert ® , Italy). A standard solution of NaOH 0.5 mol dm −3 , and the potassium dichromate, the ammonium iron (II) sulphate, and the sulfuric acid (96% v/v) used for the titrations were from Sigma-Aldrich.

Instrumentation
The multi-element analysis was performed on a NexION 300X ICP-MS spectrometer, equipped with a S10 autosampler, a glass concentric nebulizer, a glass cyclonic spray chamber, and a kinetic energy discrimination (KED) collision cell, all produced by Perkin Elmer (Milan, Italy). Samples were mineralized by means of a microwave single reaction chamber (SCR) system (ultraWAVE™, Milestone, Sorisole, Italy) equipped with fifteen polytetrafluoroethylene (PTFE) vessels (volume: 15 cm 3 each). An Ultraturrax mixer model T18 (IKA, Staufen, Germany) was used to homogenize the honey samples before analysis. The residual acidity and the dissolved organic carbon (DOC) in the digested sample were determined by acid-base titration procedure and by the Walkley-Black

Reagents and Standard Solutions
In all the analytical phases, type I water (resistivity > 18 MΩ cm −1 ), produced by means of a MilliQ plus System (Millipore, Vimodrone, Italy) was used. Nitric acid (67-69% w/w, NORMATON ® for ultra-trace analysis) and hydrogen peroxide (30% w/w, NORMATON ® for ultra-trace metal analysis) were from VWR (Milan, Italy). The multi-element standard periodic table mix 1 for ICP (TraceCert ® , 33 elements, 10 mg dm  Italy). A standard solution of NaOH 0.5 mol dm −3 , and the potassium dichromate, the ammonium iron (II) sulphate, and the sulfuric acid (96% v/v) used for the titrations were from Sigma-Aldrich.

Instrumentation
The multi-element analysis was performed on a NexION 300X ICP-MS spectrometer, equipped with a S10 autosampler, a glass concentric nebulizer, a glass cyclonic spray chamber, and a kinetic energy discrimination (KED) collision cell, all produced by Perkin Elmer (Milan, Italy). Samples were mineralized by means of a microwave single reaction chamber (SCR) system (ultraWAVE™, Milestone, Sorisole, Italy) equipped with fifteen polytetrafluoroethylene (PTFE) vessels (volume: 15 cm 3 each). An Ultraturrax mixer model T18 (IKA, Staufen, Germany) was used to homogenize the honey samples before analysis. The residual acidity and the dissolved organic carbon (DOC) in the digested sample were determined by acid-base titration procedure and by the Walkley-Black method [65] using a Thermo Scientific Orion 950 titrator, whereas the total organic carbon (TOC) was determined by means of a CHN analyzer Leco 628. Nylon filters (pore diameter: 0.22 µm), polypropylene (PP) metal-free tubes, and polyethylene (PE) flasks were from VWR (Milan, Italy).

ICP-MS Method Assessment, Quality Control and Assurance
The minimization of the dissolved organic carbon (DOC) in digested honey samples is an essential prerequisite to make a reliable ICP-MS analysis. Failing in this, high residual amounts of organic substances from a partial oxidation of the saccharides contained in honey could cause the formation of molecular ions in the plasma that may interfere with the determination of very important elements such as 52 Cr, 63 Cu, and 75 As. Beyond this specific issue, the need to avoid the interference of molecular ions of any origin on the determination of the chosen isotopes of the analytes suggested to operate in kinetic energy discrimination (KED) mode for 15 elements out of 23. Therefore, the He flow has been carefully optimized for each element, to find the best compromise between the minimization of the polyatomic ion interferences and the maximization of the instrumental signal [66][67][68][69]. As a result, 75 As, 138 Ba, 111 Cd, 59 Co, 52 Cr, 63 Cu, 57 Fe, 55 Mn, 60 Ni, 121 Sb, 120 Sn, 88 Sr, 120 Te, 51 V, and 66 Zn were analyzed in KED mode using a He flow rate between 3 and 4 cm 3 min −1 , while 107 Ag, 9 Be, 209 Bi, 7 Li, 202 Hg, 98 Mo, 208 Pb, and 205 Tl were analyzed in normal mode. The optimized instrumental parameters and the elemental settings used for each analyte are reported in Table S3.
The matrix effect was determined comparing the slopes of the calibration function obtained in the absence (2% HNO 3 v/v aqueous solutions) or presence (digested honey samples spiked with known amounts of analyte) of the matrix [70]. The data obtained accounted for a substantial absence of any matrix effect for all the analytes considered in the whole calibration range, hence quantification has been accomplished by means of external calibration. The samples, blank-diluted when necessary to lead analyte concentrations within the relevant calibration range, were analyzed in duplicate. Data obtained as a result of a triplicate ICP-MS measurement, were blank-corrected. To compensate for any signal instability, a solution of Rh (10 µg dm −3 ) was used as an internal standard, while the reliability of the measurements was ensured by analyzing one blank and a standard solution containing each analyte (50 µg dm −3 ) every 10 samples. Memory effects between consecutive samples were eliminated by interposing a washing cycle of 60 s with a 2% HNO 3 v/v aqueous solution.

Optimization of the Composition of the Acidic/Oxidizing Mixture
After some preliminary evaluations, the digestion program (Table 1; pression, temperature and time) and water volume of the oxidizing mixture (4 cm 3 ) were kept constant, while the factors considered in the design were the sample amount (0.5-1.0 g), and the ratio between the volumes of HNO 3 and H 2 O 2 (0.5/3 cm 3 -2/1.5 cm 3 ). As a result, a two-level 2 2 full factorial design was applied to improve the composition of the oxidizing mixture and minimize the dissolved organic carbon in digested sample.
To estimate the experimental error and validate the regression model, a duplicate of the central point was added to the factorial design. The responses of the design were the residual acidity in the digested sample and the efficiency of organic matter decomposition [71] (EOMD%), which was calculated with Equation (1): where TOC is the total organic carbon, in mg kg −1 ; DOC is the dissolved organic carbon, in mg kg −1 .
The experimental matrix obtained with the two-level full factorial design is reported in Table 7 as well as the experimental plan. The results obtained with the full factorial design are reported in Table 7. Multilinear regression (MLR) provided the coefficients for Equation (2) for both the responses, residual acidity (Y 1 ) and EOMD% (Y 2 ), respectively, which are reported in Table 8.
where b 0 is a constant, b 1 and b 2 are the coefficients of the main effects of the factors X 1 and X 2 , whereas b 12 is the coefficient of their interaction. As reported in Table 8, all the coefficients were significant for both the responses, including the coefficient of the interactions (b 12 ) especially for EOMD%. From the experiments at the central point, the experimental error can be estimated, and the standard deviation for Y 1 and Y 2 were 0.06 and 0.42, respectively. The predicted values for the residual acidity and for EOMD% at the center point were 0.5 ± 0.1 and 96 ± 1. Since the difference values between the experimental and predicted values are not significant (t-test, α = 0.05), the model can be accepted.
The effects of the sample amount and the ratio between the oxidizing compounds on the residual acidity and the EOMD% of the digested samples are shown in the contour plots in Figure 4.  From Figure 4 it is possible to understand the interaction between the two variables: Y1. Residual Acidity: the effect of the ratio HNO3/H2O2 (X2) on the residual acidity is present only at the lower sample amount (X1) where its increase leads to higher acidity, while at higher amount var. X2 has no effect; on the other hand, at higher ratios (X2), an increase in the sample amount (X1) led to a higher decrease in the response, greater than the decrease occurring at lower ratios (X2). Y2. EOMD%: the effect of the ratio HNO3/H2O2 (X2) on the EOMD is present only at the higher sample amount (X1) where its increase led to a higher response, while at the lower sample amount var. X2 had no effect; on the contrary, at lower ratios (X2), an increase in the sample amount (X1) led to a decrease in the response, while at higher ratios, the sample From Figure 4 it is possible to understand the interaction between the two variables: Y 1 . Residual Acidity: the effect of the ratio HNO 3 /H 2 O 2 (X 2 ) on the residual acidity is present only at the lower sample amount (X 1 ) where its increase leads to higher acidity, while at higher amount var. X 2 has no effect; on the other hand, at higher ratios (X 2 ), an increase in the sample amount (X 1 ) led to a higher decrease in the response, greater than the decrease occurring at lower ratios (X 2 ). Y 2 . EOMD%: the effect of the ratio HNO 3 /H 2 O 2 (X 2 ) on the EOMD is present only at the higher sample amount (X 1 ) where its increase led to a higher response, while at the lower sample amount var. X 2 had no effect; on the contrary, at lower ratios (X 2 ), an increase in the sample amount (X 1 ) led to a decrease in the response, while at higher ratios, the sample amount had no effect.
It is evident that a good compromise between the two responses provides a sample amount of approximately 0.7 g of honey and a ratio HNO 3 /H 2 O 2 less than 0.5. Therefore, in order to minimize the consumption of nitric acid for the benefit of the greenest hydrogen peroxide, the acceptable composition of the acidic/oxidizing mixture is 0.5 cm 3 HNO 3 , 3 cm 3 of H 2 O 2 and 4 cm 3 of H 2 O.

Statistical Analysis
A two-tail t-test at α = 0.05 was used in ascertaining the existence of matrix effect as well as in the trueness evaluation. Principal components analysis, LDA, and MLR were performed by means of the R-based software Chemometric Agile Tool (CAT) developed by the Italian group of Chemometrics [72].

Conclusions
For the first time, the concentration of 23 trace elements were measured in a very large sampling of the most renowned unifloral honeys from Sardinia, Italy, using an original and validated ICP-MS method. Special attention was paid to the development of the acid microwave digestion procedure as well as the optimization of instrumental parameters to improve the efficiency of the organic matter decomposition and to minimize the polyatomic interferences, respectively. Among the most abundant elements (i.e., Ba, Mn, Fe, Cu, and Zn), only Mn was measured in all of the samples, whereas the others ranged from the relevant LoDs to a few mg kg −1 . Toxic elements were almost always below the amounts of potential health concern, hence confirming a very good level of food safety for Sardinian honeys. Since no elemental signatures were reported in the literature for asphodel and thistle honeys, a meta-analysis was carried out only on eucalyptus and strawberry tree honeys and it highlighted the possibility of a geographical discrimination thanks to their elemental signature, mainly for strawberry tree honeys from Croatia and Sardinia. Finally, through the elemental signature of the four unifloral honeys here considered, a good classification based on the botanical origin was accomplished by means of linear discrimination analysis.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules27062009/s1, Figure S1: Scree plot; Figure S2: Influence plot; Figure S3: Permutation test. Table S1: Whole data set of the determination of 23 trace elements in 133 Sardinian unifloral honeys (Excel file); Table S2: Correlation coefficients; Table S3: Instrumental parameters and elemental settings used for the ICP-MS determination of 23 trace elements in unifloral honeys.