Harvest year effects on Apulian EVOOs evaluated by ¹H NMR based metabolomics

Sampling

A total of 900 extra virgin olive oils were obtained from four Apulian cultivars (Coratina, Cima di Mola, Ogliarola, Peranzana) in the provinces of Bari and Foggia (Apulia, Italy), collected in two subsequent harvesting periods: 2013/14 (H13) and 2014/15 (H14) (450 samples/year) (Table 1). The cultivars were initially recorded for each tree as assessed by each producer, and successively genetically confirmed for all samples by extraction of genomic DNA from leaves and detection of specific microsatellite markers, according to a standardized procedure (Salimonti et al., 2013). Each olive oil sample was produced from the drupes of a single tree (Piccinonna et al., 2016). Olive harvest were performed from each marked tree (with identification code) at optimal olive ripening stage in different periods, depending on cultivar and growing conditions (dry or well-watered field). Starting from Ogliarola cultivar, which has a natural early olive fruit ripening, the olive harvest operations were activated in the first decade of December and conducted up to the end of January for both the two harvest years (2013/14 and 2014/15). About 25–30 kg of olives per tree were collected by the use of suitable net catchers under the tree canopy, with the help of both mechanical and manual pickers and successively stored in airy boxes, marked with the identification code of the tree. Oil micro-extraction was normally performed within 24 h after the olive harvest. EVOO samples were obtained using the Spremoliva C30 milling machine (Toscana Enologica Mori, Tavarnelle Val di Pesa (FI), Italy), and successively stored in sealed dark glass bottles at room temperature in the dark prior to analysis. As much as 25–30 kg of olives were processed for each working cycle of approximately 2–3 h. The machine was cleaned carefully after each cycle.

NMR measurements

NMR samples were prepared dissolving ∼140 mg of olive oil in CDCl₃ and adjusting the mass ratio of olive oil:CDCl₃ to 13.5%:86.5%. Next, 600 µL of the prepared mixture were transferred into a 5-mm NMR tube. This ratio was chosen to give the best tradeoff for sensitivity/solution viscosity in spectral acquisition (Bruker Italia, standardized procedure for olive oil analysis). ¹H NMR spectra were recorded on a Bruker Avance spectrometer (Bruker, Karlsruhe, Germany), operating at 400.13 MHz, T = 300 K, equipped with a PABBI 5-mm inverse detection probe incorporating a z axis gradient coil. NMR experiments were performed under full automation for the entire process after loading individual samples on a Bruker Automatic Sample Changer (BACS-60), interfaced with the software IconNMR (Bruker). Automated tuning and matching, locking and shimming, and calibration of the 90°  hard pulse P(90°) were done for each sample using standard Bruker routines, ATMA, LOCK, TOPSHIM and PULSECAL, to optimize NMR conditions. For each sample, after a 5-min waiting period for temperature equilibration, a standard one-dimensional (¹H ZG) NMR experiments was performed. The relaxation delay (RD) and acquisition time (AQ) were set to 4 s and ∼3.98 s, respectively, resulting in a total recycle time of ∼7.98 s. FIDs were collected into time domain (TD) = 65,536 (64 k) complex data points by setting: spectral width (SW) = 20.5524 ppm (8223.685 Hz), receiver gain (RG) = 4, number of scans (NS) = 16. Accumulation of 16 scans (or even fewer) are usually used for samples where metabolites are present in high concentrations, as in the case of olive oil (Barison et al., 2010; Del Coco et al., 2016).

¹H NMR spectra pre-processing and Statistical analysis

The NMR raw data set was pre-processed using Topspin 2.1 and AMIX 3.9.15 (Bruker BioSpin GmbH, Rheinstetten, Germany). The FIDs were multiplied by an exponential line broadening function (0.3 Hz) before Fourier transformation and automatically phased. Spectra were referenced to the TMS signal at 0.00 ppm, used as an internal standard, obtaining good peak alignment. NMR spectra were processed using Topspin 2.1 (Bruker) and visually inspected using Amix 3.9.15 (Bruker, Biospin). Furthermore, spectra were segmented in rectangular buckets of fixed 0.04 ppm width and integrated, using the Bruker Amix software. Bucketing was performed within 10.00–0.5 ppm region, excluding the signal of the residual non-deuterated chloroform and its carbon satellites (7.6–6.9 ppm); total sum normalization was applied to minimize small differences due to total olive oil concentration and/or acquisition conditions among samples. The Pareto scaling method, which is performed by dividing the mean-centered data by the square root of the standard deviation, was then applied to the variables (Gallo et al., 2014; Sundekilde, Larsen & Bertram, 2013). The data table generated by all aligned buckets row reduced spectra was used for multivariate data analysis. Each bucket row represents the entire NMR spectrum, and all the molecules present in the sample. Each bucket in a buckets row reduced spectrum is labeled with the value of the central chemical shift for its specific 0.04 ppm width. The variables used as descriptors for each sample in chemometric analyses are the buckets. Multivariate analyses (MVA) and graphics were obtained using Simca-P version 14 (Umetrics, Sweden) using different procedures: PCA, PLS-DA and OPLS-DA (Lindon, Nicholson & Holmes, 2011). Principal Components Analysis (PCA), an unsupervised pattern recognition method, was performed to examine the intrinsic variation in the data set. To maximize the separation between sample classes, Partial Least-Squares Discriminant Analysis (PLS-DA), was applied. The PLS-DA is the regression extension of PCA, which gives the maximum covariance between the measured data (X variable, matrix of buckets related to metabolites in NMR spectra) and the response variable (Y variable, matrix of data related to the class membership). Beside PLS-DA, also Orthogonal Partial Least-Squares Discriminant Analysis (OPLS-DA) has been applied in MVA. As shown in several metabolomics recent studies, OPLS-DA represent the most recently used technique for the discrimination of samples with different characteristics (such as cultivars and/or geographical origin). OPLS-DA is a modification of the usual PLS-DA method which filters out variation that is not directly related to the response. The further improvements made by the OPLS-DA in MVA resides in the ability to separate the portion of the variance useful for predictive purposes from the not predictive variance (which is made orthogonal). Furthermore, OPLS-DA focuses the predictive information in one component, facilitating the interpretation of spectral data. On other hand, when a four categories (the cultivars) model was used for further classification purposes PLS-DA rather than OPLS-DA was preferred (Boccard & Rutledge, 2013). Both for PLS-DA and OPLS-DA, the quality of the models obtained was assessed by R² and Q² values. The first (R²) is a cross validation parameter defined as the portion of data variance explained by the models and indicates goodness of fit. The second (Q²) represents the portion of variance in the data predictable by the model. This latter indicates the model predictive ability, which is extracted according to the internal 7 fold cross-validation method of SIMCA-P software (Holmes et al., 2008; Trygg & Wold, 2002). The minimal number of components required can be easily defined since R² (cum) and Q² (cum) parameters display completely diverging behaviour as the model complexity increases. The addition of further unnecessary components to the model can, therefore, easily be detected and avoided.

Cluster difference and quality metrics

Principal Component Analysis (PCA) finds the principal components of data, provides a general overview and underlies the structure of the data. However, by visual inspection of PCA scoreplots only a qualitative general separation of the cluster can be evaluated. In order to quantify the magnitude of the separation of the clusters in the PCA scoreplots, the Mahalanobis distances were calculated (Mahalanobis, 1936). Mahalanobis distances calculated between groups in PCA scores space will closely approximate those calculated on the original data while avoiding possible collinearity of the original variables (Anderson et al., 2008). Mahalanobis distances account for different variances in each direction (PC1, PC2, PC3) and are scale-invariant. Moreover, the quality of the clusters in the number of selected PCs was estimated by using the J₂ criterion. This parameter provides a measure of the compactness and identity of the cluster (Worley, Halouska & Powers, 2013).

The J₂ criterion is defined as: ${\displaystyle J_2=\frac{\left|S_w+S_b\right|}{\left|S_w\right|};}$ with S_w being the within-class scatter and S_b the between-class scatter. A high J₂ value means well separated and tight clusters. The cluster distances and the cluster quality criterion J₂ were calculated with the open source (GNU General Public License 3.0). PCA-utils (Worley, Halouska & Powers, 2013) software, freely available on https://github.com/geekysuavo/pca-utils, and compiled on a Windows 10 64 bit notebook with the mingw-w64 gcc compiler (http://mingw-w64.org/doku.php) in the msys2 environment (http://msys2.github.io/). The data feed to the software were exported from the PCA analyses performed on Umetrics Simca v. 14 (Umetrics Software, Sweden).

Pairwise Mahalanobis distances

Data exported in table format from Umetrics SIMCA software were successively analyzed with the R statistical environment, Version 3.2.4, on a 64 bit Windows machine (R Development Core Team, 2008), using the RStudio environment, version 0.99.893 (RStudio Team, 2015). Pairwise Mahalanobis distances were calculated with the biotools package (Anderson et al., 2008). Plots were prepared with the ggplot2 library (Wickham, 2009). Some of the data preparation operations were performed with the dplyr package (Wickham & Francois, 2015).

Chemicals

All chemical reagents for analyses were of analytical grade. Deuterated chloroform (CDCl₃ 99.8%-d) containing tetramethylsilane TMS (0.03% v/v) was purchased from Armar Chemicals (Döttingen, Switzerland).

Influence of olive cultivar

As a first attempt, in order to reveal a general data grouping of all the samples, an unsupervised PCA analysis was applied to the whole data (¹H NMR-bucket-reduced spectra), revealing the presence of some outliers (22 out of 900 samples, divided over the two years), which have been excluded from the analyses. In the PCA analysis two components explained 72.6% of total variance (44.1%, 28.5% for t[1] and t[2] respectively) describing the samples distribution in the space. Visual inspection of t[1]/t[2] PCA scoreplot, reported in Fig. 1, showed a certain degree of separation in particular for the Coratina samples, that were found essentially at negative values of t[1] component and along t[2] component [0.25–0.26]. A relevant degree of overlap was observed among the three remaining classes, Cima Di Mola, Ogliarola and Peranzana. The exclusion of the Coratina class from the PCA model allowed further separation among the three remaining cultivars, especially for Cima di Mola and Ogliarola (see Fig. S3).

In order to improve the separation among the classes for all the four studied cultivars, supervised PLS-DA and OPLS-DA analyses were performed. In these methods, the identity of each sample group is specified in the model such that the maximum variance of the groups can be attained in the hyperspace. Two performance indicators were used to assess the supervised model complexity and eventual over fitting degree: the cross validation (CV) and the response permutation test (n = 400). The PLS-DA model resulted in six components with R²X = 0.88, R²Y = 0.625 and Q² = 0.61. The 2D t[1]/t[3] and 3D t[1]/t[2]/t[3] PLS-DA scoreplots showed, in this case, a clear separation for Coratina and Peranzana, and a partial overlapping for Cima di Mola and Ogliarola groups (Figs. 2A and 2C). By examining the loadings (Fig. 2B) of the original variables it was possible to define the molecular components distinctive for each class (cultivar). The Coratina group showed high values of monounsaturated fatty acids (i.e., oleic acid), as shown by corresponding loadings at δ_H 5.34, 2.02, 1.30, 1.98; the Peranzana class was characterized by high values of polyunsaturated fatty acids (PUFA), as indicated by the loadings at δ_H 2.74, 2.78, signals of linoleic and linolenic bis-allylic groups respectively. High relative content of saturated fatty acids (δ_H 1.26, corresponding to the methylene of the saturated acyl group) was found for the Cima di Mola and Ogliarola classes, partially overlapped in the t[1]/t[3] scoreplot.

Influence of harvest year

In order to obtain further information on the behaviour of the four cultivars within the studied timespan, unsupervised PCA (see Figs. S4 and S5) and supervised PLS-DA analysis was performed considering the two harvesting years separately.

The first PLS-DA model, obtained from ¹H NMR spectra of 2013/14 harvesting year EVOO samples, (six components with R²X = 0.88, R²Y = 0.625 and Q² = 0.61) showed a samples distribution in the t[1]/t[3] scoreplot (Fig. 3A) which is analogue to that observed in the PLS-DA scoreplot when considering all the samples (from both the two studied harvesting years) (see Fig. 2A). A clear separation of the Coratina group from the remaining classes was still observed especially along the first PLS-DA component (t[1]) (Fig. 3A). A certain degree of separation, in particular on the third component (t[3]), was also observed for Peranzana group, while Ogliarola and Cima di Mola appear considerably overlapped in the scoreplot. Analogously for the PLS-DA t[1]/t[3] scoreplot of the 2014/15 harvesting year, Coratina EVOO samples clearly separated from the other three cultivars (Fig. 3B). In this case, the three remaining classes are characterized by a higher degree of separation, along t[3] component (with respect to the 2013/14 samples). This result suggests that, for the studied cultivars, a potential harvesting year effect seems to influence essentially the metabolic profiles of EVOOs other than Coratina.

In order to evaluate the potential effect of the harvesting year on the EVOOs characteristics, we compared the two harvesting years for each cultivar performing unsupervised and supervised analysis on the ¹H NMR reduced spectra. The PCA analyses of the two years samples for each cultivar are reported in Figs. 4A, 5A, 6A and 7A. The four cultivars were then analyzed by OPLS-DA in order to increase the differences observed in the PCA analysis (Figs. 4C, 5C, 6C and 7C). In the case of the Coratina class, the unsupervised PCA analysis (Fig. 4A) did not reveal differences between the two periods. A poor separation between the two years could be also observed by performing a supervised OPLS-DA analysis, giving a model in which one predictive and six orthogonal components gave R²X = 0.90, R²Y = 0.85 and Q² = 0.82 (Fig. 4C). The predictive component explained 5.4% of the total variance and the uncorrelated (orthogonal) components to[1], to[2], to[3], to[4], to[5] and to[6] corresponded to 34.8%, 7.92%, 18.1%, 11.6%, 10.8% and 2.31% of the explained variance, respectively. Interestingly, by examining the loadings (Fig. 4D) of the original variables a higher relative content of saturated fatty acids (δ_H 1.26, corresponding to the methylene of the saturated acyl group) was observed for the 2013/14 Coratina samples.

In the case of the Cima di Mola samples, the unsupervised PCA analysis revealed a certain degree of separation (Fig. 5A) between the two harvest years. This separation was also improved by the supervised OPLS-DA analysis, that gave a good model with one predictive and five orthogonal components (with R²X = 0.92, R²Y = 0.95 and Q² = 0.93) (Fig. 5C). The predictive component explained 51.1% of the total variance and the uncorrelated (orthogonal) and components to[1], to[2], to[3], to[4], to[5] and to[6] corresponded to 11.7%, 13.3%, 9.8%, 11.6%, 3.8% and 2.2% of the explained variance, respectively. A nice partition of the samples, coming from the two different harvesting years, was clearly observed. A high content of saturated fatty acids (loadings at δ_H 1.26, corresponding to the methylene of the saturated acyl group) in the samples of the 2013/14 harvesting period was shown from the analysis of the loadings of the original variables, determining the separation between the two groups (Fig. 5D).

A clear partition of the samples was also observed in the case of Ogliarola class. As for the Ogliarola class, the visual inspection of t[1]/t[2] PCA scoreplot highlighted the separation among the two groups of samples (Fig. 6A). Also in this case, the OPLS-DA analysis gave a good model, (one predictive and six orthogonal components, R²X = 0.91, R²Y = 0.97 and Q² =0.95) (Fig. 6C), improving the separation between the classes. The predictive component explained 21.5% of the total variance and the uncorrelated (orthogonal) and components to[1], to[2], to[3], to[4], to[5] and to[6] corresponded to 27.4%, 19.7%, 15.7%, 3.3%, 1.5% and 1.6% of the explained variance, respectively. Interestingly, a reversal trend for the content of saturated fatty acids was observed for the Ogliarola with respect to Coratina and Cima di Mola samples, showing in this case a relatively higher content of saturated fatty acids (δ_H 1.26, corresponding to the methylene of the saturated acyl group) for the 2014/15 with respect to 2013/14 samples. Moreover, a high content of linolenic acid (δ_H 0.9, 2.06, 2.78, 5.38) in the samples of the 2013/14 harvest period was observed from the loadings of the original variables.

In the case of the Peranzana group, the unsupervised PCA analysis revealed a separation between the two harvest periods (Fig. 7A) and this separation was highlighted by the OPLS-DA analysis (Fig. 7C). The good OPLS-DA model resulted in 1 predictive and six orthogonal components, that give R²X = 0.89, R²Y = 0.94 and Q² = 0.91. The predictive component explained 24.8% of the total variance and the uncorrelated (orthogonal) components to[1], to[2], to[3], to[4], to[5] and to[6] corresponded to 31.6%, 15.4% , 6.6%, 6%, 2.8% and 2.6% of the explained variance, respectively. The separation between the two harvesting periods was due, as also observed for Ogliarola, to the higher content of linolenic acid (δ_H 0.9, 2.06, 2.78, 5.38) in the 2013/14 samples. On the other hand, similarly to Coratina and Cima di Mola samples, a higher content of saturated fatty acids (δ_H 1.26, corresponding to the methylene of the saturated acyl group) was observed by examining the S-line loading plot for the 2014/15 with respect to 2013/14 samples (Fig. 7D).

Fatty acid profile in EVOOs is primarily determined by its cultivar, although environmental factors can influence oil quality. Climatic conditions, seasonal weather fluctuations, such as rainfall and temperature can affect the physiological behavior of olive tree and the metabolic profile of its oil (Lombardo et al., 2008; Romero et al., 2003). However the effect of the seasonal variability are strongly cultivar dependent as observed in previous investigations (Inglese et al., 2011; Romero et al., 2003). We observed that Ogliarola samples showed a higher relative content of linolenic acid in 2013/14 (warm year) with respect to 2014/15 harvest year, the last characterized by dramatically intense rainfall and cool temperatures (Consorzio Lamma Rete Toscana, 2016). At the same time, Peranzana group was observed to exhibit a similar high relative content of linolenic acid. In the case of the Cima di Mola samples, a strong content of saturated fatty acids was observed in the 2013/14 harvest year. Harvest year and in particular climatic conditions resulted a discriminating parameter although Coratina EVOOs appear to be the less affected and more stable. From this point of view, the Coratina-based EVOO appears to provide well-defined chemical and sensory characteristics of a genetic and territorial origin (soil and climate) (Fanizzi et al., 2015). The lower harvesting year effect which affects the Coratina with respect to the other cultivars considered in the present study is also evident when considering the Q² predictivity parameter for the OPLS-DA models of Figs. 4B, 5B, 6B and 7B. Indeed, due to the partial overlap of 2013/2014 and 2014/2015 samples, the prediction ability for the year classification related to Coratina (Q² = 0.82) is about 10% lower with respect to Ogliarola (Q² = 0.95), Cima di Mola (Q² = 0.93) and Peranzana (Q² = 0.91) cultivars. The results obtained with this semi-quantitative assessment on the influence exerted by the harvesting year on the studied cultivars was been further examined by two additional metrics, namely the Mahalanobis distances and the J2 criterion, calculated for the four cultivars using the data from the PCA analyses shown in Figs. 1, 4A, 5A, 6A and 7A.

Cluster difference and quality metrics: Mahalanobis distances and J₂ criterion

The graphical representation of the unsupervised PCA, while giving a good relative idea of group dispersion and distances, cannot be used to directly compare among independently conducted analyses. Moreover, when the amount of cumulative variance needed to show the difference in groups requires more than three PCs, an intuitive graphical representation is not possible. Therefore, we calculated different metrics on clusters in order to obtain quantitative values useful for summarizing our conclusions. The first metric we calculated to give a numerical reference value for the PCAs was the Mahalanobis distance (MAH) among the groups. To the best of our knowledge, few chemometric studies discussed the most appropriate method on estimating the suitable number of components (PCs) to measure this type of distance (Brereton & Lloyd, 2016). For this reason, the Mahalanobis distances were measured according to two different principles: either by keeping the explained variance constant or by selecting a fixed number of components. First of all, Mahalanobis distances were calculated between the two years for the same cultivar in the PCA including all the cultivars (Table 2 with reference to the PCA of Fig. 1), and then between the two years for each single cultivar (Table 3 with reference to the PCAs of Figs. 4A, 5A, 6A and 7A).

A general overview can be obtained from the MAH distances calculated for all the PCA models studied and reported in Tables 2 and 3. In particular, when the PCA containing all the cultivars (99% cumulative variance and 20 PCs) was used, the Coratina and Ogliarola cultivars showed the relatively smallest and greatest variation, respectively, for the EVOO obtained in the two successive harvesting campaigns (Table 2, column 1). The same trend and comparable MAH values were obtained in the case of the single cultivar PCA models where, in order to obtain the 99% explained variance, a different number of components for each cultivar was considered (Table 3, column 1). Interestingly, the MAH distances calculated within the two harvesting years and with a single general PCA (Table 2, column 1) or four different cultivar specific PCAs (Table 3, column 1) all fall in the same order of magnitude, with the highest observed value (Ogliarola) less three times that of the lowest (Coratina). It should be noted that the distances calculated only with the first 2 PCs, both for a single general PCA (Table 2, column 2) or four different cultivar specific PCAs (Table 3, column 2) could be misleading, since the amount of variance explained by the first 2 PCs is not close, therefore making the distances not comparable. In fact, Mahalanobis distances calculated for the single PCA first two PCs (Table 2, column 2) showed Ogliarola and Peranzana having a comparable behaviour, with Cima di Mola appearing to be the most sensitive cultivar to the harvesting year related oil variation. Cima di Mola appeared to be the most sensitive cultivar to the harvesting year effect also in the case of the four different, cultivar specific, PCAs (Table 3, column 2), while Ogliarola and Peranzana showed a similar MAH. Nevertheless, even using only the first 2 PCs both for a single general PCA (Table 2, column 2) or four different cultivar specific PCAs (Table 3, column 2) the calculated MAH distances clearly show that the Coratina EVOO could be considered the less sensitive to the harvesting year related variations.

Finally, in order to measure the compactness and identity of each cluster we applied the J₂ criterion. The J₂ values were calculated for all the cultivars for both the two harvesting years, in the single PCA model built for the whole dataset. The J₂ values for each cultivar could describe also the cultivar specific internal variability in each harvesting year. A high value of J₂ indicates well separated and tightly clustered groups of samples. Although J₂ values could be also calculated for 99% variance (20 PCs) (see Table S1), for clarity reasons we have compared in Table 4 the obtained J₂ values in the cases of two or three considered PCs. In the first case, the J₂ values reported the scatter of each cluster in first two dimensions, for both the harvesting years, describing shape and relative positions for the ellipses containing each group in the bidimensional plane. Similarly, the J₂ values reported for 3 PCs described the scatter of each cluster in the first three dimensions, defining shape and relative positions for the ellipsoids containing each group in the Cartesian space. Differently from the previous metric, the J₂ values maintain a consistent behaviour, as represented in the examined cases (Table 4). The cluster integrity and separation are very similar both in the two cases (2 and 3 PCs) and for 99% variances, where the J₂ criterion is able to give a measure of the quality of clustering also in a n-dimensional space (20 PCs). Cima di Mola showed the highest J₂ values for both the two harvests, indicating well-separated compact clusters, especially for 2014/15. Ogliarola and Peranzana showed similar behaviors for both the two harvest years, with a very high J₂ value (108.67) obtained for Peranzana in the 2014/15 harvest year. Finally, low J₂ values reported for Coratina indicated a very high homogeneity of EVOO samples for both the two harvests and number of considered PCs. This data confirmed the Coratina cultivar to be the less sensitive to the harvesting year related variations, as already discussed for the MAH distances.

The absolute values increased in the case of higher dimensionality of the data, but the trends are maintained. Except for the Coratina class, the J₂ values are consistently lower for the 2013/14 with respect to 2014/15 harvest year, indicating both higher overlap and diffusion of the clusters. In Fig. 1 and Fig. S4 the J₂ values can be visually interpreted: Coratina has the widest variance (apparently due to the presence of two possible subgroups), while the Cima di Mola has the tightest cluster. Peranzana has few overlaps with the other cultivars, but at the same time has a high variance. The same results are visible qualitatively for the 2014/15 harvest, but all the clusters appear to be more spaced among themselves (Fig. S5). When the number of PCs examined goes to four, the general trends are maintained, with the spread for the cultivars increasing, except for the Coratina.

Pairwise Mahalanobis distances

As largely described in literature (Jonsson et al., 2015; Westerhuis et al., 2010), NMR-based metabolomics is a very powerful techniques but affected by different types of bias, for example when a large intrinsic variation is present between the subjects or data acquired over long periods of time were combined and analyzed. This is particularly evident in clinic and diagnostic studies, when intrinsic variation between human and/or biological samples can largely influence metabolome. For this reason, multivariate analysis, and in particular the supervised OPLS-DA technique, works reducing the uncorrelated and/or unwanted variation by applying a separation between the subjects variation from the studied effect, focusing the attention on the effect of interest. In our work, we applied this criterion when considering the whole dataset and the clustering quality of unsupervised models. In the next and last step of this work, a paired (rather than random) structure of the data was considered (Fig. 8), in order to provide information about the difference across the study population.

To the best of our knowledge, there are no other reports on harvesting year effect for EVOO samples originated from individual olive trees with genetically established cultivar identification.

Using the PCA data, the pairwise squared Mahalanobis distances (MADs) were calculated for all the couples of samples (2013/14 and 2014/15 harvest year). The samples having no correspondence in the other harvest (due to harvesting or oil production) were excluded from this analysis. The square roots of said values were calculated and then divided in two sets: a set consisting of the distances between the EVOO obtained in the two harvest years from a single plant, and a set containing the distances of a single sample in the harvest 2013/14 from all the samples of the same cultivar in the following harvest (2014/15). The modulus of the distances was considered. The first distribution can be envisioned as a trace of the difference in the global environment that a single plant experienced in the two harvest year, while the second distribution could be held to represent an all-including representation of the cultivar variability between the two harvests.

Notwithstanding the different sample sizes for the four cultivars, an information that can be extracted by the data is that the couples of cultivars Coratina–Ogliarola and Cima di Mola-Peranzana share a similar behavior. In fact, Coratina and Ogliarola have a smaller mean distance and median between the two harvests (Fig. 9 and Table 5), indicating that on average the samples had a lower variability in the two harvests than Cima di Mola and Peranzana, although Coratina has a much higher number of measured samples, so a higher variability would not have been unexpected. On the other hand, standard deviations are close for all the cultivars except a slightly higher value for Coratina, as expected on the basis of the much wider sample size. Cima di Mola and Peranzana share higher medians and means, with a more symmetric distribution around the mean. These data suggest a higher variability of the EVOO obtained from a single plant in the two years, a possible consequence of significant differences in the climatic sub regions where the plant were located, between two harvest years.

The trends in the aggregated data are similar (Fig. 10 and Table 6). The overall shape of the distributions and the average distances remain close, with Peranzana featuring slightly smaller mean and median. This could be attributed to the change in spread of the Peranzana in Harvest 2014/15, with the formation of two subgroups, which can be easily visually identified in Figs. 7A and 7B, confirmed by the doubling of the J₂ value of H14 with respect to H13.