Combination of Stable Isotope Analysis and Chemometrics to Discriminate Geoclimatically and Temporally the Virgin Olive Oils from Three Mediterranean Countries

The knowledge of the isotopic composition of virgin olive oil (VOO) allows the evaluation of authenticity and geographical origin, being an important tool against fraud. This study aimed to assess if VOOs produced in three Mediterranean regions could be discriminated on the basis of multivariate statistical analysis of geoclimatic and isotopic data. A total of 138 geo-referenced VOO samples from Portugal, France and Turkey from two different cultivation years were collected. The isotopic composition (δ13C, δ2H and δ18O) of VOOs was obtained using an elemental analyzer coupled to an isotope ratio mass spectrometer (EA-IRMS). One-way analysis of variance for δ13C, δ2H and δ18O showed some significant differences either between crop years or geoclimatic conditions. Based on multiple regression analyses using meteorological and geographical parameters, a meteoric water line for olive oil from Portugal, France and Turkey, in two harvest years, were created to assess the impact of climate change on their δ2H and δ18O values. Principal component analysis and Linear Discriminant Analysis, used to sort samples according to geoclimatic origin, performed best for French and Portuguese olive oils. In light of the results, multivariate isotopic analysis of VOO samples may discriminate not only between geoclimatic regions but also among cultivation years.


Introduction
Traceability of olive oils is still a popular research topic, addressed in many different ways, including exploring the benefits of different analytical tools due to the increased demand for an effective analytical methodology that can verify the geographic origin of olive oils. As claim of origin is becoming more important, so is the need for new tools to establish it. Traditional chromatography techniques (gas or liquid) coupled to mass spectrometry [1][2][3] or UV detection [4], or the "more environmentally friendly" analytical methods based on spectroscopy [5,6] have been used for this purpose.
However, olive oil stable isotope ratio analysis (SIRA), either in bulk, using Isotope Ratio Mass Spectrometry (IRMS) coupled to an elemental analyzer (EA), or in a class of compounds using IRMS Foods 2020, 9,1855 2 of 12 coupled with a gas chromatographer (GC), is nowadays used more often to achieve separation among olive oils from different geographical origins [7][8][9][10][11][12]. The SIRA is based on the fact that the stable isotope content of bio-elements reflects the geo-climatic characteristics of the area of production [13] and was initially performed in bulk analysis. In the case of olive oil, bulk δ 13 C, δ 18 O and δ 2 H determinations allowed the differentiation among different geographic origins, at least in some cases. With innovations in this analytical technique that now allow the study of stale isotopes in specific compounds, like fatty acids, or volatile compounds, more detailed information can be retrieved allowing improvement in geographical discrimination.
Oxygen and hydrogen isotope ratios (δ 18 O and δ 2 H), as a natural tracer, have been effectively used in diagnosing water sources [14], as well as a marker of authenticity and origin of food or beverages. Particularly, these isotopes have been employed to trace the authenticity and geographical origin of European extra-virgin olive oils [12,15]. The stable isotope compositions of oxygen and hydrogen in soil and ground waters are controlled by local precipitation in the recharge and catchment areas, depending on various geoclimatic factors (e.g., temperature, distance from the sea, amount of atmospheric precipitation, etc.) [16,17]. Therefore, the variations in δ 2 H and δ 18 O values in plants and their products, through their relationship with humidity-precipitation and temperature conditions during plant growth, may allow the tracing of the geographic origin of vegetable oils [18][19][20][21].
Bearing this in mind, our study aimed to assess if VOOs from three different Mediterranean countries (Portugal, France and Turkey) could be discriminated on the basis of multivariate statistical analysis of geoclimatic and isotopic data (δ 13 C, δ 2 H and δ 18 O) obtained from bulk analysis of 138 varietal VOOs. For the majority of the selected VOO samples, approximately 5 kg of each different olive variety was processed separately in an Abencor ® system (Seville, Spain) within 24 h after harvesting. Fruits were crushed with a hammer mill and the olive paste was malaxed at 25 • C for 30 min in an olive paste mixer. Finally, the olive oil was separated by centrifugation. Other samples were taken directly from olive oil mills to avoid possible undeclared mixtures with olive oils from other cultivars and geographical origins, before bottling. They were stored in dark-brown glass bottles at 20 • C in the dark. The VOO samples were perfectly geo-referenced, obtaining data for latitude (UTM), longitude (UTM), altitude (m.a.s.l.), sea distance (km), mean annual temperature ( • C) and mean annual rainfall (mm) (Table S1).

Stable Isotope Analysis
The carbon, hydrogen and oxygen isotope composition ratios ( 13 C/ 12 C, 2 H/ 1 H and 18 O/ 16 O, respectively) of VOO samples were determined by elemental analysis/isotope ratio mass spectrometry (EA/IRMS) following the procedure described by Jiménez-Morillo et al. [12]. Briefly, bulk isotopic signature of light elements (δ 13 C, δ 2 H and δ 18 O) was analysed using a Flash 2000 HT with two reactors (combustion and pyrolysis, respectively). The elemental analyzer is coupled by a ConFlo IV interface unit to a continuous flow Delta V Advantage isotope ratio mass spectrometer (IRMS) (Thermo Scientific, Bremen, Germany). Isotopic ratios "δ" are reported in milli Urey (mUr) from appropriate standards recognized by the International Atomic Energy Agency (IAEA). The standard deviation of bulk δ 13 C, δ 2 H and δ 18 O were ± 0.1, 1.0 and 0.5 mUr, respectively. Each sample was measured in duplicate (n = 2).

Statistical Analysis
Multivariate statistical analysis of data was carried out with Statgraphics Centurion XV software (Statgraphics Technologies, Inc, The Plains, VA, USA). Principal component analysis (PCA) was used for simultaneous ordination of different geographic and climatic variables and the δ 13 C, δ 2 H and δ 18 O values, illustrating their mutual relationships. Multiple linear regression (MLR) was applied to generate an isotopic forecasting model by harvest year, using geographic and climatic factors as independent variables. The validation of the model was confirmed by the production of spurious models computed from the fully randomized dependent variables (isotopic composition of hydrogen, oxygen and carbon) from the 138 VOO samples from three different Mediterranean countries in two harvest years. The alternative models displayed a poor correlation (p > 0.05) with geoclimatic variables (data not shown). Therefore, the well-correlated models generated with experimental data were not overfitted [12,22]. Linear discriminant analysis was also carried out to evaluate whether olive oil samples from the three countries could be easily separated based on isotopic data. The statistical significance of each discriminant function was evaluated by Wilk's lambda [23] and the predictive ability of the LDA model was evaluated by leave-one-out cross validation.

Stable Isotope Analysis of VOO Samples
Table S1 collects the carbon, hydrogen and oxygen isotope values of VOO samples from production areas in France, Portugal and Turkey in two harvest years. Concerning δ 13 C value, the Portuguese VOO samples ranged between −31. These isotopic values are in accordance with those reported by Gumus et al. [9] for olive oils from western Turkey. The δ 13 C values of the VOO samples from the three Mediterranean countries were within the typical value of C 3 photosystem plants [24]. However, as observed, there were significant differences between location and harvest years. One-way analysis of variance (one-way ANOVA) using δ 13 C composition ( Figure 1a) as independent variable (IV) indicated that Portuguese and Turkish VOOs are significantly (p < 0.05) different between themselves, but not between respective cultivation years. The French samples were significantly different (p < 0.05) between cultivation years. Traditionally, it has been accepted that δ 13 C composition of food-stuff depends on the localization of the production area, i.e., on the geomorphological and climatic conditions, such as temperature, rainfall, sea distance, longitude and latitude [12]. These factors influence leaf stomata opening and hence the photosynthesis yield [25], as dry and warm conditions may produce a restriction in the opening time of stomas, limiting the admission of atmospheric CO 2 to the leaf, causing an increase in δ 13 C value. The existence of no significant differences between harvest year in Portugal and Turkey may be related to the similitude in mean annual temperature between 2016 and 2017 (Table S2). However, in French regions a significant difference between harvest year was observed, which was also displayed in the temperature values (Table S2).
On the other hand, hydrogen and oxygen isotope composition of VOO samples (Table S1) (Table S1), the values ranged: (i) between −146.6 and −127.1 mUr, and between −155.6 and −122.0 mUr for Portugal; (ii) between −186.7 and −141.8 mUr, and between −178.2 and −127.1 mUr for France; and (iii) between −159.0 and −125.6 mUr, and between −178.1 and −135.1 mUr for Turkey. One-way ANOVA analysis of δ 2 H values ( Figure 1c) showed that only the samples from Turkey were significantly (p < 0.05) different according to the cultivation year. In contrast, the samples from Portugal and France are significantly (p < 0.05) different from each other, which may be directly related to a quite different climate, mainly temperature and precipitation (Table S2). With respect to Turkey, the significant differences between harvest years may be due to specific climatological events that could have altered the isotopic composition of hydrogen. This hypothesis may be supported by the fact that the δ 2 H and δ 18 O composition of foodstuffs are significantly correlated to (1) the local climatic conditions (temperature and precipitation) [12,16,26]; (2) soil and plant evapotranspiration [27,28]; (3) air humidity; and (4) fractionation associated with the biosynthesis of plant tissues [29,30]. Therefore, the differences observed between countries and harvest years may be related to the above mentioned factors. However, the existence of irrigation activities may also modify the isotopic composition of hydrogen and oxygen [31,32].  Figure 2 shows the correlation between δ 2 H and δ 18 O of VOO samples from the three Mediterranean countries in two harvest years. In agreement with other researchers, the H and O isotope values of VOO samples are directly correlated to the meteoric water uptake by olive plants, so it has been observed that the H isotope value of olive oil is about 100 times lower than the value of rainwater, while the O isotope value is about 27 times higher than that of rainwater [30,33]. Therefore, in this work, these criteria were used to generate the local meteoric water line (LMWL) for each location and harvest year ( Figure 2). In addition, the correlation lines were compared with the global meteoric water line (GMWL). The GMWL was generated using the δ 2 H and δ 18 O values obtained from IAEA WISER (Water Isotope System for data analysis visualization and Electronic Retrieval database) [34]. The slope variation of each LMWL with the GMWL allows assessment of the alteration of climatic conditions in a geographical location and year [35,36]. In general, the LMWL for the Mediterranean countries studied displayed a depletion of slope in comparison with the GMWL, which is related to rise in temperature and reduction of precipitation [35]. The current climatic change and global warming scenario is dominated by the rising of mean temperature in the Mediterranean basin together with the reduction of the number of rainfall episodes and the amount of water [37,38]. A more detailed analysis of δ 2 H and δ 18 O by country and year shows that in 2016 all studied countries showed the same depletion slope, which may be correlated to a rise in drought conditions. However, in 2017, the trends were different. In the case of French samples in 2017, a slope parallel to the GMWL is displayed, while in 2016 its slope displays a conspicuous depletion. The average value of temperature and precipitation (Table S2) in France showed a significant (p < 0.05) depletion of 39% in rainfall in 2017. A possible explanation may be the existence of an irrigation process, since, in this year, the cultivation area suffered a significant reduction in rainfall (over 30%, Table S2). Irrigation may alter the natural Foods 2020, 9, 1855 5 of 12 δ 2 H and δ 18 O value [39,40]. Therefore, the non-alteration of the 2017 slope in France may be due to the anthropogenic factor. Displacement may be seen on the y-axis (δ 18 O composition), which may be related to the isotope irrigation water value. In the case of Portugal, slopes in different harvest years display the same value. However, in 2017, the LMWL showed a shift in δ 18 O composition, i.e., in the y-intercept. This may be supported by a significant reduction in rainfall of 31% (Table S2). On the other hand, Turkey shows a slight change in slope in 2017. The analysis of precipitation data (Table  S2) showed a rise in rainfall of 21%, but this is not significant (p > 0.05). The good correlation of LMWL (correlation between treated δ 2 H and δ 18 O of VOOs), as well as their connection to the GMWL, conditioned by inherent geoclimatic factors (temperature, latitude, sea distance, rainfall, etc.), confirm the potential use of isotope composition of VOOs for their geoclimatic characterization. Some recent works have explored this possibility in VOO samples produced from the same olive variety and harvest year in different countries, for example, in Italy [19] and Spain [20].

Meteoric Water Line
Foods 2020, 9, x FOR PEER REVIEW 5 of 12 depletion of 39% in rainfall in 2017. A possible explanation may be the existence of an irrigation process, since, in this year, the cultivation area suffered a significant reduction in rainfall (over 30%, Table S2). Irrigation may alter the natural δ 2 H and δ 18 O value [39,40]. Therefore, the non-alteration of the 2017 slope in France may be due to the anthropogenic factor. Displacement may be seen on the y-axis (δ 18 O composition), which may be related to the isotope irrigation water value. In the case of Portugal, slopes in different harvest years display the same value. However, in 2017, the LMWL showed a shift in δ 18 O composition, i.e., in the y-intercept. This may be supported by a significant reduction in rainfall of 31% (Table S2). On the other hand, Turkey shows a slight change in slope in 2017. The analysis of precipitation data (Table S2) showed a rise in rainfall of 21%, but this is not significant (p > 0.05). The good correlation of LMWL (correlation between treated δ 2 H and δ 18 O of VOOs), as well as their connection to the GMWL, conditioned by inherent geoclimatic factors (temperature, latitude, sea distance, rainfall, etc.), confirm the potential use of isotope composition of VOOs for their geoclimatic characterization. Some recent works have explored this possibility in VOO samples produced from the same olive variety and harvest year in different countries, for example, in Italy [19] and Spain [20].

Chemometric Analysis
Principal component analysis (PCA) was employed on the one hand to identify the possible correlation between independent (stable isotopes of VOO samples) and dependent (geoclimatic parameters of produced areas) variables (Figure 3a), and on the other hand to sort the VOO samples by geographical region and harvest year (Figure 3b). The analysis of the three first components (component 1: 29.35%, component 2: 20.95%, and component 3: 17.37%) can explain up to 67% of the total variance (Table S3). The scatterplot of the loadings of PC-1 vs. PC-2 vs. PC-3 (Figure 3a) showed that δ 2 H values were directly correlated with temperature, sea distance and altitude, and negatively with latitude, longitude and rainfall. In contrast, δ 18 O values were correlated with longitude, altitude and temperature, while longitude and rainfall were negatively correlated with oxygen isotope composition. These correlations were previously observed in other articles by the authors. For

Chemometric Analysis
Principal component analysis (PCA) was employed on the one hand to identify the possible correlation between independent (stable isotopes of VOO samples) and dependent (geoclimatic parameters of produced areas) variables (Figure 3a), and on the other hand to sort the VOO samples by geographical region and harvest year (Figure 3b). The analysis of the three first components (component 1: 29.35%, component 2: 20.95%, and component 3: 17.37%) can explain up to 67% of the total variance (Table S3). The scatterplot of the loadings of PC-1 vs. PC-2 vs. PC-3 (Figure 3a) showed that δ 2 H values were directly correlated with temperature, sea distance and altitude, and negatively with latitude, longitude and rainfall. In contrast, δ 18 O values were correlated with longitude, altitude and temperature, while longitude and rainfall were negatively correlated with oxygen isotope composition. These correlations were previously observed in other articles by the authors. For example, Jiménez-Morillo et al., [11], using for the first time the analytical pyrolysis-compound specific isotope analysis (Py-CSIA), found that temperature and sea distance were directly correlated with the hydrogen isotope composition of specific compounds of EVOO samples from five Mediterranean countries. However, other authors showed inversely correlation between δ 18 O and sea distance [12]. With respect to carbon isotope (δ 13 C), this displays direct correlation with longitude and altitude but correlates inversely with rainfall. Camin et al. [41] and Jiménez-Morillo et al. [12] observed the same correlations between carbon isotope composition of olive oil samples and the geoclimatic factors. More rain will cause plant's stomata to stay open longer, causing more isotopically light CO 2 molecules to pass inside the leaf and the δ 13 C value to drop, relative to plants growing in a more arid environment. The existing correlations observed in Figure 3a were corroborated by a multi linear regression (MLR) approach. The MLR, using the geoclimatic parameters of the VOO producer locals (n = 138) from the studied Mediterranean countries, led to significant (p < 0.05) forecasting models for the carbon, hydrogen and oxygen isotope composition for the harvest year.
Foods 2020, 9, x FOR PEER REVIEW 6 of 12 example, Jiménez-Morillo et al., [11], using for the first time the analytical pyrolysis-compound specific isotope analysis (Py-CSIA), found that temperature and sea distance were directly correlated with the hydrogen isotope composition of specific compounds of EVOO samples from five Mediterranean countries. However, other authors showed inversely correlation between δ 18 O and sea distance [12]. With respect to carbon isotope (δ 13 C), this displays direct correlation with longitude and altitude but correlates inversely with rainfall. Camin et al. [41] and Jiménez-Morillo et al. [12] observed the same correlations between carbon isotope composition of olive oil samples and the geoclimatic factors. More rain will cause plant's stomata to stay open longer, causing more isotopically light CO2 molecules to pass inside the leaf and the δ 13 C value to drop, relative to plants growing in a more arid environment. The existing correlations observed in Figure 3a were corroborated by a multi linear regression (MLR) approach. The MLR, using the geoclimatic parameters of the VOO producer locals (n = 138) from the studied Mediterranean countries, led to significant (p < 0.05) forecasting models for the carbon, hydrogen and oxygen isotope composition for the harvest year.    This cross-validation method has been previously used in several works [11,12,42]. The MLR analysis allowed us to obtain the Equations (1)-(6) of the prediction model, using automatic backwards variable selection: In this way, the new predicting isotope value models (Equations (1)-(6)) will allow forecasting of changes that foodstuffs will undergo in different climate change scenarios, and adapting the criteria for the protection of designations of origin, safeguarding society from possible food fraud. In addition, they could be used to assess the isotopic composition of virgin olive oils under different climatic change scenarios (from the mildest to the most severe). However, it is known that the isotopic forecasting models developed should be corrected every year owing to intrinsic fluctuations in isotopic composition. Therefore, the next step will be the validation and correction of the possible isotopic alterations, generating correction factors.
On the other hand, the 3D scores in the PCA plot of VOO samples (Figure 3b), generated using their unique isotopic composition, displayed three clusters, one for each country. In the case of the French and Turkish cluster, their separation is dominated by the difference in the stable isotope value of hydrogen and oxygen. The French samples displayed higher isotope values than the VOO samples from Turkey (Figure 2b,c). The Portuguese cluster showed an intermediate location in the 3D PCA plot (Figure 3b). These samples show no significant difference with respect to the other two regional VOO samples (Figure 2b,c). Although the samples are well-grouped into three clusters depending on the country, a more detailed study shows that the samples can be separated slightly depending on the harvest year. The most prominent case is the samples from Portugal (Figure 3b). This difference may be directly related to the variation in the δ 18 O (Figure 2b), as well as the significant difference in rainfall (Table S2). With respect to the other Mediterranean VOO samples, this difference between cultivation years was not so clear.
In order to verify if a better classification of olive oils can be achieved, a Linear Discriminant Analysis (LDA) was also performed. Figure 5 shows the scatter plot for olive oils from the three countries in each year. Two discriminant functions were constructed based on Wilks' Lambda values, which explained 100% of the variance, 81.2% of the total variance explained by function 1 and 18.8% explained by function 2, for 2016 olive oils; and 95.3% of the total variance explained by function 1 and 4.7% explained by function 2, for 2017 olive oils (Table S5). The Wilks' Lambda values were, for function 1 and 2, 0.377 and 0.795 with a p-value of 0.000 and 0.001, for 2016 olive oils; and 0.230 and 0.878 with a p value of 0.000 and 0.010, for 2017 olive oils, which implies that there are significant differences among the three regions and the first function has higher discrimination power, since Wilk's Lambda is lower and p-value is equal to 0. Discriminant functions can be shown as:  Table S5 summarizes the performance of the discriminant models regarding the classification of olive oil samples, and shows that the application of cross-validation (CV) exposes some misclassifications, mainly regarding samples from Turkey. For 2016 samples, 79.0% of the olive oil was correctly classified and this value diminished to 74.2% after cross validation; for 2017 samples, results were better, with 92.9% of the cases correctly classified before and after CV. Olive oils from Portugal performed better than the others with the highest percentage of correctly classified olive oils, according to country of origin.

Conclusions
In light of the results, the combination of multivariate isotopic analysis of VOO samples and chemometrics allows discrimination between geoclimatic regions but also among cultivation years. Therefore, it presents a powerful tool to fight against food fraud in the VOO sector. In addition, the design of this work has allowed the generation of isotopic forecasting models using the geoclimatic parameters as predictive factors. Therefore, the interpolation of geoclimatic data (mainly temperature and rainfall) under different climatic change scenarios (from benevolent to drastic) will allow us to assess the isotopic alteration of virgin olive oil in future scenarios.
This displays a great potential for establishing new correction factors to determine the geographical origin (Designation of Origin) of olive oils.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Table S1: Varieties, origin, region, geo-climatic information and stable isotope composition of VOO samples; Table S2: Average and Standard deviation of Temperature and Rainfall in France, Portugal and Turkey in two years (2016 and 2017). Alteration of rainfall amount (%); Table S3: Principal component loading scores; Table S4: Variable standardized canonical discriminant function coefficients along with discriminant functions classification parameters; Table  S5: Overall classification performance (count and %) of the LDA model, with and without cross validation (CV), for the three countries differentiation study.  Table S5 summarizes the performance of the discriminant models regarding the classification of olive oil samples, and shows that the application of cross-validation (CV) exposes some misclassifications, mainly regarding samples from Turkey. For 2016 samples, 79.0% of the olive oil was correctly classified and this value diminished to 74.2% after cross validation; for 2017 samples, results were better, with 92.9% of the cases correctly classified before and after CV. Olive oils from Portugal performed better than the others with the highest percentage of correctly classified olive oils, according to country of origin.

Conclusions
In light of the results, the combination of multivariate isotopic analysis of VOO samples and chemometrics allows discrimination between geoclimatic regions but also among cultivation years. Therefore, it presents a powerful tool to fight against food fraud in the VOO sector. In addition, the design of this work has allowed the generation of isotopic forecasting models using the geoclimatic parameters as predictive factors. Therefore, the interpolation of geoclimatic data (mainly temperature and rainfall) under different climatic change scenarios (from benevolent to drastic) will allow us to assess the isotopic alteration of virgin olive oil in future scenarios.
This displays a great potential for establishing new correction factors to determine the geographical origin (Designation of Origin) of olive oils.
Supplementary Materials: The following are available online at http://www.mdpi.com/2304-8158/9/12/1855/s1, Table S1: Varieties, origin, region, geo-climatic information and stable isotope composition of VOO samples; Table S2: Average and Standard deviation of Temperature and Rainfall in France, Portugal and Turkey in two years (2016 and 2017). Alteration of rainfall amount (%); Table S3: Principal component loading scores; Table S4: Variable standardized canonical discriminant function coefficients along with discriminant functions classification parameters; Table S5: Overall classification performance (count and %) of the LDA model, with and without cross validation (CV), for the three countries differentiation study.