Rare Earth Element Variability in Italian Extra Virgin Olive Oils from Abruzzo Region

Extra virgin olive oil is a food product from the Mediterranean area that is particularly and continuously experiencing to increasing instances of fraudulent geographical labeling. Therefore, origin protection must be improved, mainly based on its intrinsic chemical composition. This study aimed to perform a preliminary chemical characterization of Abruzzo extra virgin olive oils (EVOOs) using rare earth elements (REEs). REEs were evaluated in EVOO samples of different varieties produced in different geographical origins within the Abruzzo region (Italy) in three harvest years using ICP-MS chemometric techniques. Principal component, discriminant, and hierarchical cluster analyses were conducted to verify the influence of the variety, origin, and vintage of the REE composition. The results of a three-year study showed a uniform REE pattern and a strong correlation in most EVOOs, in particular for Y, La, Ce, and Nd. However, europium and erbium were also found in some oil samples. Compared with cultivar and origin, only the harvest year slightly influenced the REE composition, highlighting the interactions of the olive system with the climate and soil chemistry that could affect the multielement composition of EVOOs.


Introduction
Extra virgin olive oil (EVOO) is an important agricultural product.Its consumption is common in everyday cooking [1,2].It represents the main fat source of the "Mediterranean diet", providing a special aroma and taste due to the contribution of its phytochemical compounds [3].Its high contents of monounsaturated fatty acid and antioxidant compounds has beneficial health effects [4][5][6] in preventing coronary diseases, cancer types, diabetes, and autoimmune illness [7][8][9].Worldwide, olive oil production was estimated as almost 3 million tons for the period from 2021 to 2022, with about 2 million tons produced in Europe (Spain, 1400 t; Italy, 329 t; Greece, 232 t; and Portugal, 206 t) and about 1 million tons in non-European countries (Tunisia, 240 t; Turkey, 235 t; Morocco, 200 t; Algeria, 91 t; Egypt, 20 t; and Argentina, 3 t) [10].
The Abruzzo region produces the fifth largest EVOO volume in Italy, at about 144,000 tons in 2021 [11].Owing to local varieties and the specific pedoclimatic conditions in the Abruzzo region, it is possible to cultivate olive trees over the entire territory, starting from the sea to the foothills of Majella and the Gran Sasso mountains located 600-700 m above sea level [12].Fifty percent of the annual regional EVOO production is concentrated in the Chieti district, and the Pescara province produces thirty percent.The rest of the production is distributed in Teramo and L'Aquila districts, with 16 and 4%, respectively.The EVOO produced in the country is appreciated all over the world as a very-high-quality product [13], and, for this reason, tracing the origin of EVOO to control product authenticity or adulteration is fundamental [14][15][16][17][18].In fact, the adulteration of EVOO with other vegetable oils is relatively easy to recognize.Conversely, it is difficult to check extra virgin oils produced from foreign countries or those using imported olives and sold as local products, because the overall composition of the final product can be almost totally similar to that of local ones [19].
Here, the focus is not only on food safety and quality control but also on the assessment of the authenticity of the declared geographical origin.Indeed, frauds can have a significant impact on customer health, confidence in the product, and, consequently, e final consumption behavior, causing significant economic losses [20,21].For these reasons, EC Regulation 1151/2012 provides rules on marketing standards for olive oil and states the mandatory nature of origin labeling [22].It introduces discrimination based on geographical indicators, highlighting that the characteristics of the product are also related to the geographical origin [23].The Protected Designation of Origin (PDO) label on olive oil has become a motivating choice criterion for customers, since olive oil quality and flavor are linked to the origin of the olives, and they are associated with specific production practices [24,25].
Even though a method for certifying EVOOs' geographical origin has not yet been established in the scientific literature, the discrimination of olive oil production on a geographic basis, also called geographical authentication, is becoming an important trend in scientific research, with different analytical chemistry approaches being applied, including elemental/isotopic, nuclear magnetic resonance (NMR), mass spectroscopy (MS), and energy-dispersive X-ray fluorescence (XRF), as well as organic analytical methods and organoleptic evaluation [22,[26][27][28].The organic analytical methods are based on the determination of fatty acids and triacylglycerols [29], but minor metabolites such as phenols [30], aliphatic and terpene alcohols, sesquiterpene hydrocarbons [31], sterols [32,33], and pigments are also considered.The isotopic profile analysis method is mostly based on the detection of stable isotope (H, C, and O) ratios, since the isotopic fractionation is correlated with geographical and climatic parameters [33,34].The isotopic composition of strontium (Sr) is also useful, since it reflects local soil and geology [35][36][37],and, in particular, the geogenic soil formation that overlies the geological substrate [38,39].
Alternatively, authentication can be realized by evaluating the mineral composition of olive oils, which gives information regarding the biological demand of the plant, the bioavailability and mobility of mineral compounds from the soil, and the influence of agronomic practices such as the use of fertilizers and pesticides.Since olive trees of the same cultivar, characterized by the same genetics but planted in different countries, produce different oils, the analysis of their elemental profiles reflects the effect of their interaction with pedoclimatic conditions and agricultural practices [40][41][42].For this reason, to achieve precise authentication, it is necessary to identify elements acting as natural markers that are not influenced by secondary sources such as farming practices and anthropic sources and are suitable for the determination of the origin of EVOO [26,43].
Trace elements concentration in the Earth's crust do not exceed 1 g/kg, while, in olive oil, their concentrations, except for Ca, do not exceed a few hundred micrograms per kilogram (µg/kg) [27].The concentrations in olive oil may differ according to several parameters such as olive cultivar [44], irrigation of olive trees [45], use of fertilizers [46], fluctuation in annual climatic parameters [47], and the positive fractionation that occurs during the active absorption from the soil and the translocation to the fruits [19].For these reasons, no consensus has been reached with regard to a direct and clear correlation between their concentration in the soil and in the related olive oils [19,[47][48][49].
Foods 2024, 13, 141 3 of 16 Ultra-trace elements, also called rare earth elements (REEs), show wide soil concentration variability around the world, ranging from 0.1 to 700 mg kg −1 [50,51].Conversely, several authors have found REE concentrations in olive oils as ranging from 0.002 to 7 ng g −1 [52].Ultra-trace elements are not essential for the growth and development of olive tree and are passively absorbed without active fractionation.Therefore, the REE composition can provide a representative fingerprint of the olive oil since REEs more proportionally reflect their less-abundant concentrations in soils.Conversely, fractionation from the original distribution in soil occurs when a certain element is assumed in a preferential way because it is a nutrient, or it is excluded from absorption because it is toxic [19,26,36,41,[53][54][55][56][57].
For these reasons, REEs are suitable for discriminating foodstuffs on a geographic basis, acting as geochemical markers as well [58], while major, minor, and trace elements are useful in authentication schemes based on varietal or technological discrimination [19,26].Joebstl and coauthors used REEs to identify the geographical origin of pumpkin seed oil [59] and, lately, some studies deepened the understanding of the relationship among the production chain and EVOO characteristics [19,57].Barbera et al. confirmed that the relationship between the soil and olive fruits depends exclusively on the soil REE composition [56] and identified an excellent marker to identify the geographical fingerprint of EVOOs [57].
Based on the previous assumptions, the aim of this study was to evaluate, for the first time, the REE contents of EVOOs from the Abruzzo region, mainly from Chieti and Pescara provinces, which produce 80% of the total EVOO in the region.In particular, 29 EVOOs produced during three harvest years (2019, 2020, and 2021) were fingerprinted with inductively coupled plasma mass spectrometry (ICP-MS).The dataset was analyzed using chemometric tools, in particular, using principal component analysis (PCA), linear discriminant analysis (LDA), and hierarchical clustering analysis (HCA), with the aim of discriminating between olive oil groups.

Area Sampled, Local Geology, and Sample Collection
Olive harvest and EVOO extraction were conducted during the 2019, 2020, and 2021 harvesting years.The locations of the olive orchards where the drupes were collected are reported in Figure 1, which shows a simplified geological map of the Abruzzo region and its position in the coastal Adriatic region (i.e., the outer (eastern) portion of the Apennines orogenic belt (geological map of Italy: 1:50,000 scale, 361 sheets).In this region, carbonate rocks in the mountain range are exposed and sandstones alternate with marls in the adjacent low-lying mountain-foot domains.The soils of the Adriatic coastal region have developed above a common sedimentary substrate made by clayey-sandy-conglomerate shallow marine deposits of the Mutignano Formation (Late Pliocene-Lower Pleistocene).Locally, the soil is featured by silty-clayey eluvial-colluvial, alluvial, and slope Quaternary continental deposits.All the investigated olive orchards, therefore, were located on soils that had pedogenized the areas of the Mutignano Formation (FMT), with the exception of the Alanno site, which, however, developed on lithologies very similar to those of the FMT.
exception of the Alanno site, which, however, developed on lithologies very similar to those of the FMT.Olives were mechanically harvested in October from nonexperimental orchards, from 20 plants for each cultivar.Olive fruits were at a medium level of ripening, which is specific for each variety and defined by each producer.EVOOs were extracted within 12-24 h of harvest in an olive mill located in the same municipality as the olive orchard, with industrial, three-phase, continuous extraction systems.Oils were collected during the extraction, packaged in dark 750 mL glass bottles, and stored in dark conditions at about 18 °C until analysis.Oil samples were coded with alphanumeric codes based on the harvest year (Table 1); "A", "B", and "C" corresponded to 2019, 2020, and 2021, respectively.The trees and industrial olive mills were the same during the three years of the research.Information about the geographical location of the orchards, sample coding, and cultivars are presented in Table 1.Olives were mechanically harvested in October from nonexperimental orchards, from 20 plants for each cultivar.Olive fruits were at a medium level of ripening, which is specific for each variety and defined by each producer.EVOOs were extracted within 12-24 h of harvest in an olive mill located in the same municipality as the olive orchard, with industrial, three-phase, continuous extraction systems.Oils were collected during the extraction, packaged in dark 750 mL glass bottles, and stored in dark conditions at about 18 • C until analysis.Oil samples were coded with alphanumeric codes based on the harvest year (Table 1); "A", "B", and "C" corresponded to 2019, 2020, and 2021, respectively.The trees and industrial olive mills were the same during the three years of the research.Information about the geographical location of the orchards, sample coding, and cultivars are presented in Table 1.

Olive Oil Preparation
Prior to inductively coupled plasma mass spectrometry (ICP-MS) analysis, the organic matter contained in the olive oils was destroyed via mineralization.The EVOO mineralization was carried out based on the EN 13805:2014 protocol adopted by the Istituto Zooprofilattico Sperimentale dell'Abruzzo e del Molise "G.Caporale" (Teramo, Italy).Briefly, a microwave digestion system, Milestone ultraWAVE ECR (Sorisole, Italy), equipped with temperature and pressure control was used to digest samples.About 300 mg of each EVOO sample was weighed directly into disposable glass vessels.The vessels were then filled with 4 mL of HNO 3 (60%) for trace analysis (Merck KGaA, Darmstadt, Germany) and prepared for digestion.Each digestion cycle was programmed according to the protocol: step 1-ramp to temperature 230 • C and pressure 150 bar in 25 min, at 1500 W power; step 2-the same conditions maintained for 10 min; step 3-cooling cycle.Each sample resulting from acid digestion was then diluted to 15 mL with high-purity water (18.2MΩ cm −1 resistivity) obtained from an ELGA LabWater PURELAB Option-Q water purification system (High Wycombe, UK).

Abruzzo Province Geographical Location (Number of Samples) Sample Code Cultivar
Prefixes "A", "B", and "C" in the sample code correspond to 2019, 2020, and 2021, respectively.Numbers of samples are reported in brackets (n).
The instrument used was an Agilent 7900 ICP-MS (Agilent Technologies, Tokyo, Japan), which was used in the Laboratory of Newborn Screening, Proteomics and Endocrinology of CAST, University of Chieti.Detailed operating conditions and instrumental parameters are given in Table S1.The 4th-generation Octopole Reaction System (ORS) was used to measure, at the same time, some elements (Na, Mg, K, Ca, Fe, Zn, Rb, Al, Ga, Cr, Pb) in helium (He) mode to reduce spectral interference and noise effects; other elements (V, Mn, Sr, Ba, and REE) were measured in no-gas mode due to the lack of interference [63].
The optimization of ICP-MS was carried out to obtain the maximum signal intensities for 7Li, 89Y, 140Ce, and 205Tl using a 1 µg L −1 tuning solution containing Li, Y, Co, Ce, Mg, and Tl (Agilent Technologies, Palo Alto, CA, USA), while keeping the formation of oxides 140 CeO + /140Ce + and doubly charged species Ce 2+ /Ce + ratios below 1% and 2%, respectively.The sample-introduction system was washed between analyses with 2% HNO 3 .Two multielement mixtures at 10 µg mL −1 were used in acid solution: (A) Ag, Ba, Be, Cd, Co, Cr, Cu, Fe, Mn, Ni, Pb, Rb, Se, Sr, Tl, U, V, and Zn in 5% HNO 3 ; (B) Ce, Dy, Er, Eu, Ga, Gd, Ho, In, La, Lu, Nb, Nd, Pr, Sm, Th, Tb, Tm, Y, and Yb in 5% HNO 3 .These mixtures were employed to prepare diluted calibration solutions daily, and three calibration curves were prepared using these multielement mixtures.An internal standard correction was performed via the online addition of an internal standard solution of Rh (20 µg L −1 ) in a T piece.Ultra-pure water (18 MΩ cm −1 ) was obtained from a Milli-Q system (Millipore, Bedford, MA, USA).Nitric acid (69% v/v) and an internal standard solution of Rh were bought from Merk (Darmstadt, Germany) and were ultrapure-grade.Duplicate analysis was performed for each sample.The full data were recorded with Agilent MassHunter Data Acquisition software (version 4.2) and processed with Agilent MassHunter Data Analysis software (version 4.2).
The limit of detection (LOD) and limit of quantification (LOQ) were calculated for the REEs by analyzing ten experimental blanks.The LOD and LOQ were initially calculated as signals by employing Equations ( 1) and ( 2), respectively: where t is the constant from a one-sided Student's t-test at the 95% confidence level for n − 1 degrees of freedom, yb is the average blank signal, and sb is the corresponding standard deviation.The corresponding LOD and LOQ concentration values (Table S2) were obtained by using an appropriate calibration curve satisfying the following relationship: 0.5x 1 < LOD < x 1 , where x 1 is the concentration of the first calibration level [64].Data below the LOD and LOQ values were excluded from the analysis.

Data Analysis
Empirical analyses were conducted to study the characteristics of the REEs according to some variables of interest, and the main features of the observed data were summarized by performing an overall descriptive analysis.The differences between mean REE concentrations related to year, variety, and origin were evaluated via analysis of variance (ANOVA) and Kruskal-Wallis tests.Where the ANOVA test assumptions were not satisfied (normality and homoscedasticity of residues), the Kruskal-Wallis nonparametric test was used.Significant differences were established at the level of p = 0.05.
In the second step, commonly used foodomics chemometric techniques such as PCA and HCA were used as exploratory methods, without any a priori knowledge of groups present in the population [22].To reduce collinearity among data, PCA was applied to REE concentrations, and, based on the first two principal components, the LDA, an a priori knowledge of group membership technique, was performed [65].Finally, the HCA method, revealing groups of similarity (clusters), was conducted using Euclidean distances and Ward's linkage methods.The graphical output of the analysis was a heatmap, a tree-like plot, where both rows and columns were clustered [66,67].

Elemental Profile of Olive Oils
The mean REE concentrations detected in the EVOOs are shown in Figure 2a, while the mean REE contents (≥LOQ) of each olive oil sample during the 2019-2021 period from different geographical origins are reported in Table 2. Despite the limited data on REE contents available in the literature, the concentration range that we detected ranged between 0.09 and 4.22 ng g −1 , in accordance with the results of different authors, which were well summarized in a recent study [69], who reported ranges varying between 0.002 and 7 ng g −1 in olive oils from costal region.Our values also agree with those described by other authors in other areas [2,[70][71][72].In descending order, the most abundant REEs were Ce, Y, La, and Nd, which had mean values of 3.27, 2.18, 1.53, and 1.36 ng g −1 , respectively.The levels of the remaining elements, such as Pr, Sm, Gd, Dy, Er, Yb, and Eu ranged from 0.38 to 0.13 ng g −1 .Table S3 provides a detailed list of the REE concentrations found in the literature.To better study and visualize the REE patterns, chondrite-normalized values were used [73], which provide a reference for the normalization of rare earth elements, since they are assumed to reflect the original composition of the Earth's crust [74,75].Except for Eu and Er, which presented slight positive and negative anomalies, respectively, all the oil samples showed very similar chondrite normalized patterns, confirming the above-mentioned hypothesis regarding the limited fractionation that occurs when passing from soil to fruits and to the final EVO [19].It is evident that the homogeneous lithological and geomorphological context (Figure 1) strongly affected the distribution patterns of the REEs in soils and, therefore, the general uniform REE pattern observed in the EVOOs.Slight content variations, moreover, could reflect the local influence of bedrock-weathering and soil-leaching processes.
The behavior of REE distributions strictly followed the Oddo-Harkins rule, which states how elements with an even atomic numbers are more abundant than elements with immediately adjacent atomic numbers [76,77].The REE concentration patterns of the EVOOs (Figure 2b) therefore proportionally reflect the distribution of REEs in soils, pointing out that different vintages, varieties, and locations seem did not affect the patterns.This aspect represents one of the prerequisites for a reliable chemical marker of geographical origin, and so it would be useful to further investigate the potential of REEs in this sense.

Elemental Profile of Olive Oils
The mean REE concentrations detected in the EVOOs are shown in Figure 2a, while the mean REE contents (≥LOQ) of each olive oil sample during the 2019-2021 period from different geographical origins are reported in Table 2. Despite the limited data on REE contents available in the literature, the concentration range that we detected ranged between 0.09 and 4.22 ng g −1 , in accordance with the results of different authors, which were well summarized in a recent study [69], who reported ranges varying between 0.002 and 7 ng g −1 in olive oils from costal region.Our values also agree with those described by other authors in other areas [2,[70][71][72].In descending order, the most abundant REEs were Ce, Y, La, and Nd, which had mean values of 3.27, 2.18, 1.53, and 1.36 ng g −1 , respectively.The levels of the remaining elements, such as Pr, Sm, Gd, Dy, Er, Yb, and Eu ranged from 0.38 to 0.13 ng g −1 .Table S3 provides a detailed list of the REE concentrations found in the literature.To better study and visualize the REE patterns, chondrite-normalized values were used [73], which provide a reference for the normalization of rare earth elements, since they are assumed to reflect the original composition of the Earth's crust [74,75].Except for Eu and Er, which presented slight positive and negative anomalies, respectively, all the oil samples showed very similar chondrite normalized patterns, confirming the above-mentioned hypothesis regarding the limited fractionation that occurs when passing from soil to fruits and to the final EVO [19].It is evident that the homogeneous lithological and geomorphological context (Figure 1) strongly affected the distribution patterns of the REEs in soils and, therefore, the general uniform REE pattern observed in the EVOOs.Slight content variations, moreover, could reflect the local influence of bedrock-weathering and soil-leaching processes.
The behavior of REE distributions strictly followed the Oddo-Harkins rule, which states how elements with an even atomic numbers are more abundant than elements with immediately adjacent atomic numbers [76,77].The REE concentration patterns of the EVOOs (Figure 2b) therefore proportionally reflect the distribution of REEs in soils, pointing out that different vintages, varieties, and locations seem did not affect the patterns.This aspect represents one of the prerequisites for a reliable chemical marker of geographical origin, and so it would be useful to further investigate the potential of REEs in this sense.To evaluate the strength of the relationships among the REEs, Pearson correlation coefficients were calculated and are reported in a correlation matrix (Table 3).A strong correlation (p < 0.001) was found among Y, La, Ce, Pr, Nd, Gd, and Dy; five of them were grouped in the light fraction (LREEs).Conversely, the weakly correlated REEs were Eu-Er (−0.04) and Eu-Yb (0.03).It was evident that similar chondrite patterns corresponded to stronger concentration correlations.

ANOVA and Kruskal-Wallis Test
The one-way ANOVA was conducted to evaluate the differences between mean REE concentrations related to the three vintages (2019, 2020, and 2021).The parametric test assumptions were verified for all REEs, except for Eu (failed both normality and Levene's test) and Er (failed normality test).Ten out of eleven REEs did not exhibit any significant differences (p > 0.05) since there was no variation in the mean REE concentrations among the vintages.Eu exhibited significant differences between the 2019 and 2020 vintages (according on Tukey's HSD test).
The nonparametric Kruskal-Wallis test was conducted to evaluate the differences between median REE concentrations related to the 12 varieties or cultivars (CVs) and 6 origins.No variation in the REE concentrations was found regardless of origin and the (p > 0.05), confirming the homogenous pattern distribution that we previously observed.

PCA
According to the values reported in Table 3, to reduce the collinearity among the data, principal component analysis (PCA) of the REE concentrations was conducted, and the results are shown in Table 4.The first two principal components (F1 and F2) explained 75.77% of the total variance (specifically, F1 explained 65.37% and F2 explained 10.4%).The main variables that influenced the first component were the concentrations of Y, La, Ce, and Nd (F1 eigenvectors > 0.35), while the variables Eu and Er mainly affected the second component, with eigenvector values of 0.76 and −0.54, respectively.
According to the loading plot (Figure S1), the vector lengths of Y, La, Ce, and Nd represented the largest contribution of the variance in the first component.They were located close to each other, with a small angle in between; out toward the same periphery, they covariation was strongly positive (very similar patterns) and proportional to the degree distance from the PC origin.On the other side, europium (Eu), especially, and erbium (Er), mainly contributed to the second component.The large angle between Eu and Er highlighted the weak correlation.Sm and Yb, with the shortest vector and being closest to the origin, presented the lowest absolute variances, and might be better explained by other factors [57].The first component represented elements with strong associations, reflecting similar behavior or bioavailability in the soil, belonging to the same geochemical groups on the periodic table, being able to form clusters [47,49,78].Consequently, according to Goldshmidt's geochemical classification, La, Ce, and Nd (LREEs), which have more basic behavior and higher solubility, were grouped together [47].In addition to this grouping of elements, yttrium (HREE), which also characterizes the elemental profile of olive oils, was shaped by the geochemical processes in the study territory.
A PCA biplot (Figure 3) highlights that samples C4 and C5 presented the highest concentrations of REEs and were the most positively correlated with F1.Conversely, EVOOs A4, B5, and B8 presented the lowest concentrations of REEs and were negatively correlated with F1.With respect to F2, the A samples, belonging to the 2019 vintage, specifically A4, A9, and A5, presented a highly positive correlation with europium and a negative one with erbium.Finally, the B and C samples, representing the 2020 and 2021 vintages, respectively, were mostly grouped together.According to the loading plot (Figure S1), the vector lengths of Y, La, Ce, and Nd represented the largest contribution of the variance in the first component.They were located close to each other, with a small angle in between; out toward the same periphery, they covariation was strongly positive (very similar patterns) and proportional to the degree distance from the PC origin.On the other side, europium (Eu), especially, and erbium (Er), mainly contributed to the second component.The large angle between Eu and Er highlighted the weak correlation.Sm and Yb, with the shortest vector and being closest to the origin, presented the lowest absolute variances, and might be better explained by other factors [57].
The first component represented elements with strong associations, reflecting similar behavior or bioavailability in the soil, belonging to the same geochemical groups on the periodic table, being able to form clusters [47,49,78].Consequently, according to Goldshmidt's geochemical classification, La, Ce, and Nd (LREEs), which have more basic behavior and higher solubility, were grouped together [47].In addition to this grouping of elements, yttrium (HREE), which also characterizes the elemental profile of olive oils, was shaped by the geochemical processes in the study territory.
A PCA biplot (Figure 3) highlights that samples C4 and C5 presented the highest concentrations of REEs and were the most positively correlated with F1.Conversely, EVOOs A4, B5, and B8 presented the lowest concentrations of REEs and were negatively correlated with F1.With respect to F2, the A samples, belonging to the 2019 vintage, specifically A4, A9, and A5, presented a highly positive correlation with europium and a negative one with erbium.Finally, the B and C samples, representing the 2020 and 2021 vintages, respectively, were mostly grouped together.Based on the PCA assumptions, a preliminary comparison analysis was performed including the REE data of three different Italian monocultivar oils: Pisciottana (Calabria region, D1), Frantoio (Piemonte region, D2), and Taggiasca (Liguria region, D3), as well as the data from one EU/extra EU EVOO blend (D4).The PCA score plot (Figure S2) shows the main separation of the samples along F1.EVOOs D1 and D3 were found to be similar to C4 and C5.EVOO D2 was comparable to A2 and opposite from D1 and D3.Sample D4 was positioned close to the origin.This result highlighted the better characterization of Abruzzo EVOOs and Liguria, Calabria, and Piemonte oils compared with the EVOO blend.

Discriminant Analysis
Linear discriminant analysis (LDA) was performed as a further unsupervised data elaboration, according to year, variety, and origin.It was applied to the first two principal component factors (F1 and F2), and results are depicted in confusion matrices (Tables 5, S4 and S5, Supplementary Materials).The mean correct classification rates were 78.6% for vintage class (Table 5), 46.4% for geographical location (Table S4), and 32.1% (Table S5) for the cultivar class.The vintage discrimination supported the previous PCA evaluations, identifying the A samples' distribution in the F1 direction, which were opposite to the compact B grouping, in the observation plot (Figure 4).Furthermore, samples C4 and C5 were spatially separated as evident C outliers along the F2 direction.
Based on the PCA assumptions, a preliminary comparison analysis was performed including the REE data of three different Italian monocultivar oils: Pisciottana (Calabria region, D1), Frantoio (Piemonte region, D2), and Taggiasca (Liguria region, D3), as well as the data from one EU/extra EU EVOO blend (D4).The PCA score plot (Figure S2) shows the main separation of the samples along F1.EVOOs D1 and D3 were found to be similar to C4 and C5.EVOO D2 was comparable to A2 and opposite from D1 and D3.Sample D4 was positioned close to the origin.This result highlighted the better characterization of Abruzzo EVOOs and Liguria, Calabria, and Piemonte oils compared with the EVOO blend.

Discriminant Analysis
Linear discriminant analysis (LDA) performed as a further unsupervised data elaboration, according to year, variety, and origin.It was applied to the first two principal component factors (F1 and F2), and results are depicted in confusion matrices (Table 5 and Tables S4 and S5, Supplementary Material).The mean correct classification rates were 78.6% for vintage class (Table 5), 46.4% for geographical location (Table S4), and 32.1% (Table S5) for te cultivar class.
The vintage discrimination supported the previous PCA evaluations, identifying the A samples' distribution in the F1 direction, which were opposite to the compact B grouping, in the observation plot (Figure 4).Furthermore, samples C4 and C5 were spatially separated as evident C outliers along the F2 direction.According to Table 5, the 2020 vintage samples were more accurately classified (91.7%) than those from 2019 and 2021 (75% and 62.5%, respectively).Indeed, from a graphical point of view (Figure 4), they were located closer to the reference centroid (B) compared with the A and C samples, where the distance between samples and their respective centroids was higher.
It is important to note that the low rate of LDA results related to the geographical location and the cultivar agree with the pairwise Kruskal-Wallis comparison outcomes, which did not identify significant differences among the EVOOs.

Cluster Analysis
Aggregative hierarchical cluster analysis (HCA), using Euclidean distances and Ward's linkage method, was implemented to interpret the chemometric data based on an input matrix consisting of 11 chemical variables (REEs) and 28 oil samples.The results of HCA are shown in a heatmap plot (Figure 5).

Conclusions
The overall results of this preliminary study show that in the high-EVOO-producing region of Abruzzo, the REE concentration patterns of olive oil, among different varieties, origins, and vintages, were almost homogeneous, showing the potential to be used as a marker of geographical origin.Among these three factors, only vintage slightly influenced REE concentrations, suggesting the possible effect of interactions among soil geochemisand climatic characteristics.Some REEs are more effective and useful than others in representing the transfer of soil geochemistry to olive oil, in particular, Y, La, Ce, Nd, because of their larger contribution to the overall variance, having the strongest correlations, and the most similar patterns.The research outcomes add useful data to the scarce literature on the REE concentrations in olive oils, which are important for assessing food quality, especially in Mediterranean countries, where EVOO represents the main fat source, which is consumed daily.The use of the information gathered in this study can be a useful starting point for producing a reliable discrimination model that allows the verification and certification of the geographical origin of EVOOs produced in the Chieti and Pescara districts of Abruzzo Region.Future use of the developed procedure with larger data sets will verify its value.Furthermore, the relationships between EVOO samples originating from different areas and the effects of biogeochemical drivers on the geographical distribution of the elements require further exploration.
upplementary Materials: The following supporting information can be downloaded at: www.mdpi.com/xxx/s1, Figure S1: PCA loading plot; Figure S2.Score plot of EVOO samples and foreign oils from Pisciottana, (Calabria region, D1), Frantoio (Piemonte region, D2), Taggiasca (Liguria region, D3), and EVOO blend (D4); Table S1.ICP-MS instrumentation and operating condi- The column clustering highlights the correspondence of the graphical heatmap output information with a PCA loading plot (Figure S1).Eu, Er, Sm and Yb showed independent behavior and did not group with the two main clusters of La/Ce/Nd/Y and Gd/Dy/Pr.
The row clustering analysis identified four main groups, also highlighted in the score plots (Figures 3 and 4), where C4 and C5 presented the highest concentrations of REEs, while EVOOs A4, B8, and B5, were characterized by the lowest concentrations of REEs.Moreover, EVOOs A5, A6, A7, and A9 presented the highest concentrations of Eu, while B1, B2, B10, B11, and A2 had the highest concentrations of Er.
HCA also confirmed that only the vintage influenced the REE concentrations' tendency to form groupings, while geographical location and cultivar did not, as previously highlighted by the corresponding confusion matrices.Since the average annual temperature and cumulative precipitation were similar in the three considered years (data not shown), climate discrimination should be studied considering intrayearly (monthly and seasonal) climate trends, which could influence the specific phenological stages of olive tree development.Edaphic factors, such as the bioavailability of inorganic elements in soil and its chemical characteristics (pH, electrical conductivity, organic matter and inorganic carbon (CaCO 3 ), should also be considered [71].

Figure 1 .
Figure 1.(a) Simplified geological map showing the different lithologies characterizing the Abruzzo region (modified from a geological map of Italy on a 1:50,000 scale, sheets: 361; www.isprambiente.gov.it/Media/carg/note_illustrative/361_Chieti.pdf,accessed on 09 October 2023), processed using Q-GIS tools and overlaid on a 3D topographic digital elevation model of Italy with a 10 m cell size [62].Green diamonds indicate the localization of the analyzed olive orchards: LA, Loreto Aprutino; Pi, Pianella; Al, Alanno; Ca, Casoli; Sc, Scerni; Va, Vasto.(b) Location of Abruzzo region in Italy.

Figure 1 .
Figure 1.(a) Simplified geological map showing the different lithologies characterizing the Abruzzo region (modified from a geological map of Italy on a 1:50,000 scale, sheets: 361; www.isprambiente.gov.it/Media/carg/note_illustrative/361_Chieti.pdf,accessed on 09 October 2023), processed using Q-GIS tools and overlaid on a 3D topographic digital elevation model of Italy with a 10 m cell size [62].Green diamonds indicate the localization of the analyzed olive orchards: LA, Loreto Aprutino; Pi, Pianella; Al, Alanno; Ca, Casoli; Sc, Scerni; Va, Vasto.(b) Location of Abruzzo region in Italy.

Table 1 .
Olive oil samples: geographical location, sample coding, and cultivars.

Table 1 .
Olive oil samples: geographical location, sample coding, and cultivars.

Table 2 .
REE concentrations (ng g −1 ) in olive oil samples.Mean and standard deviation (sd) of n = 3 replicates.

Table 2 .
REE concentrations (ng g −1 ) in olive oil samples.Mean and standard deviation (sd) of n = 3 replicates.

Table 5 .
Confusion matrix of samples grouped for vintage.

Table 5 .
Confusion matrix of samples grouped for vintage.