antioxidants Article Impact of Phenol-Enriched Olive Oils on Serum Metabonome and Its Relationship with Cardiometabolic Parameters: A Randomized, Double-Blind, Cross-Over, Controlled Trial Marta Farràs 1,2,3, *, Jonathan Richard Swann 4 , Ian Rowland 5 , Laura Rubió 6 , Isaac Subirana 7,8 , Úrsula Catalán 9 , Maria José Motilva 6 , Rosa Solà 9 , Maria Isabel Covas 10 , Francisco Blanco-Vaca 11,12,13,† , Montserrat Fitó 14,15,† and Jordi Mayneris-Perxachs 14,16,17, *,† 1 Unit of Nutrition and Cancer, Epidemiology Research Program, Catalan Institute of Oncology (ICO), Bellvitge Biomedical Research Institute (IDIBELL), L’Hospitalet de Llobregat, 08908 Barcelona, Spain 2 Division of Metabolism, Digestion and Reproduction, Faculty of Medicine, Imperial College, London SW7 2BX, UK 3 Institut d’Investigacions Biomèdiques (IIB) Sant Pau, C/Sant Quintí 77, 08041 Barcelona, Spain 4 School of Human Development and Health, Faculty of Medicine, University of Southampton, Southampton SO17 1BJ, UK 5 Hugh Sinclair Human Nutrition Group, University of Reading, Reading RG6 6DH, UK 6 Antioxidants Research Group, Food Technology Department, AGROTECNIO-CERCA Center, University of Lleida, Avda/Alcalde Rovira Roure 191, 25198 Lleida, Spain 7 CIBER de Enfermedades Cardiovasculares (CIBERCV), C/Monforte de Lemos 3-5, 28029 Madrid, Spain 8 Cardiovascular Epidemiology Group, IMIM (Hospital del Mar Medical Research Institute), Citation: Farràs, M.; Swann, J.R.; C/Doctor Aiguader 88, 08003 Barcelona, Spain 9 Functional Nutrition, Oxidation, and Cardiovascular Diseases Group (NFOC-Salut), Facultat de Medicina i Rowland, I.; Rubió, L.; Subirana, I.; Ciències de la Salut, Hospital Universitari Sant Joan, Universitat Rovira i Virgili, C/Sant Llorenç 21, Catalán, Ú.; Motilva, M.J.; Solà, R.; 43201 Reus, Spain Covas, M.I.; Blanco-Vaca, F.; et al. 10 NuProas Handesbolag (NUPROAS HB), 13141 Nacka, Sweden Impact of Phenol-Enriched Olive 11 CIBER de Diabetes y Enfermedades Metabólicas Asociadas (CIBERDEM), C/Monforte de Lemos 3-5, Oils on Serum Metabonome and Its 28029 Madrid, Spain Relationship with Cardiometabolic 12 Servei de Bioquímica, Hospital Santa Creu i Sant Pau—IIB Sant Pau, 08025 Barcelona, Spain 13 Departament de Bioquímica i Biologia Molecular, Universitat Autònoma Barcelona, 08193 Barcelona, Spain Parameters: A Randomized, 14 CIBER Fisiopatología de la Obesidad y Nutrición (CIBERobn), Instituto de Salud Carlos III, Double-Blind, Cross-Over, Controlled Trial. Antioxidants 2022, 28029 Madrid, Spain 15 Cardiovascular Risk and Nutrition Research Group, Regicor Study Group, IMIM (Hospital del Mar Medical 11, 1964. https://doi.org/10.3390/ Research Institute), Doctor Aiguader 88, 08003 Barcelona, Spain antiox11101964 16 Department of Diabetes, Endocrinology and Nutrition, Dr. Josep Trueta University Hospital, Academic Editors: Manolis Fousteris 17190 Girona, Spain 17 Nutrition, Eumetabolism and Health Group, Girona Biomedical Research Institute (IdibGi), and Fotini Lamari 17190 Girona, Spain Received: 10 August 2022 * Correspondence: martafarras@iconcologia.net (M.F.); jmayneris@idibgi.org (J.M.-P.) Accepted: 25 September 2022 † These authors contributed equally to this work. Published: 30 September 2022 Abstract: Phenol-rich foods consumption such as virgin olive oil (VOO) has been shown to have Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in beneficial effects on cardiovascular diseases. The broader biochemical impact of VOO and phenol- published maps and institutional affil- enriched OOs remains, however, unclear. A randomized, double-blind, cross-over, controlled trial was iations. performed with thirty-three hypercholesterolemic individuals who ingested for 3-weeks (25 mL/day): (1) an OO enriched with its own olive oil phenolic compounds (PCs) (500 ppm; FOO); (2) an OO enriched with its own olive oil PCs (250 ppm) plus thyme PCs (250 ppm; FOOT); and (3) a VOO with low phenolic content (80 ppm). Serum lipid and glycemic profiles, serum 1 H-NMR spectroscopy- Copyright: © 2022 by the authors. based metabolomics, endothelial function, blood pressure, and cardiovascular risk were measured. Licensee MDPI, Basel, Switzerland. We combined OPLS-DA with machine learning modelling to identify metabolites discrimination of This article is an open access article the treatment groups. Both phenol-enriched OO interventions decreased the levels of glutamine, distributed under the terms and creatinine, creatine, dimethylamine, and histidine in comparison to VOO one. In addition, FOOT conditions of the Creative Commons decreased the plasma levels of glycine and DMSO2 compared to VOO, while FOO decreased the Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ circulating alanine concentrations but increased the plasma levels of acetone and 3-HB compared to 4.0/). VOO. Based on these findings, phenol-enriched OOs were shown to result in a favorable shift in the Antioxidants 2022, 11, 1964. https://doi.org/10.3390/antiox11101964 https://www.mdpi.com/journal/antioxidants Antioxidants 2022, 11, 1964 2 of 16 circulating metabolic phenotype, inducing a reduction in metabolites associated with cardiometabolic diseases. Keywords: functional olive oil; phenolic compounds; metabonomics; cardiovascular diseases 1. Introduction Virgin olive oil (VOO) phenolic compounds (PCs) have been reported to possess antioxidant and anti-inflammatory properties, and exert chemoprotective effects in ex- perimental studies [1,2]. Moreover, PCs found in VOO produce beneficial changes in the serum lipid profile and haemostasis, reduce blood pressure, and have anti-thrombotic and anti-inflammatory activity in humans [3–5]. Data from the EUROLIVE study demonstrated an increase of HDL cholesterol (HDL-C), and a decrease in in vivo lipid oxidative damage, in an OOPC dose-dependent manner in healthy humans [6]. Furthermore, an enhancement in HDL function with virgin olive oil (VOO) in healthy humans has been reported [7]. In addition, VOO PCs have been reported to exert protection against risk factors for coronary heart disease, particularly in individuals with high oxidative stress [1,2]. Phenol-enriched foods can display a dual action since antioxidants may revert to pro-oxidant actors [8–10]. Functional foods with complementary PCs could therefore have their beneficial effects enhanced. In this regard, enriching OO with PCs, without increasing the fat content, could be an effective method to raise the intake of PCs and increase their potential to improve health. Our group have previously demonstrated that phenol-enriched OO improved endothelial function in pre-/hypertensive subjects [11], and enhanced endothelial and HDL functions in hypercholesterolemic individuals [12,13]. In this regard, phenol-enriched OO also increased the expression of cholesterol efflux-related genes in two transcriptomic sub-studies [14,15]. Nutritional metabolomics, or nutrimetabolomics, is being increasingly employed to study molecular interactions between diet and the global metabolic system. It allows the study of metabolic responses to dietary modulations, and permits the inter-individual variation in responses to the intake of specific nutrients and diets to be established. In addition, the metabolic phenotype contains signals derived from dietary inputs, as well as information related to the metabolic activity of the intestinal microbiota and its interactions with the diet and host. We have previously described interactions between phenol-enriched OOs and the human microbiota [16]. Food enriched with complementary antioxidants, according to their structure/activity relationship, may be a possibility to achieve healthy effects avoiding thes harmful ones. The aim of this study was to measure the serum metabolic phenotypes of hypercholesterolemic patients, and characterize the biochemical modulations induced by the intake of functional OOs enriched with their own PCs or with their own PCs plus complementary ones from thyme. 2. Materials and Methods 2.1. Phenol-Enriched Olive Oil Preparation and Characteristics A VOO with a low-phenolic content (80 ppm) was used as a control condition and as a matrix of enrichment to prepare two functional OOs (FOOs; 500 ppm). The first was a functional VOO (FOO; 500 ppm) enriched with its own PCs by the addition of a phenolic extract obtained from freeze-dried olive cake. The phenolic profile of the olive cake extract was comparable to the control VOO, as they were obtained from the same olive variety (Arbequina cv) and the same olive-growing area. The second FOO was enriched with its own PCs and complemented with thyme PCs using a phenolic extract obtained from a mixture of freeze-dried olive cake and dried thyme (Thymus vulgaris) (FOOT; 500 ppm). The FOOT thus contained 50% OOPCs (some of them were naturally present in VOO and some of them came from the added olive cake extract) and 50% thyme PCs, which have been Antioxidants 2022, 11, 1964 3 of 16 shown to possess a range of anti-inflammatory and anti-oxidative properties [17]. For the wash-out period a common OO was used. The procedure to obtain the phenolic extracts and the phenol-enriched OO has been published elsewhere [18]. Briefly, the phenolic extracts were obtained using and accelerated solvent extractor as previously optimized [18,19] obtaining an average of 2 g (1.96 ± 0.16 g) of freeze-dried phenolic extract from 10 g of raw material (freeze-dried olive cake and dried thyme). Once the phenolic extracts were obtained, they were used to prepare phenol-enriched OOs in the ratio of 2.5 g of extract/100 g of VOO with low phenolic content using a dispersion method with water (2%) until complete homogenization. The phenolic identification and quantification of the control VOO and the phenol-enriched OOs were analyzed by HPLC coupled with tandem mass spectrometry (MS/MS) as previously described [18]. As previously reported, the OOs only differed in their phenolic content, the fat and micronutrient compositions were the same (Supplementary Table S1). 2.2. Study Design The VOO and HDL Functionality (VOHF) study was a randomized, double-blind, crossover, controlled trial including 33 volunteers with cholesterol levels higher than rec- ommended (total-cholesterol > 200 mg/dL) (19 men/14 women), aged 35 to 80. Exclusion criteria were body mass index (BMI) > 35 Kg/m2 , smokers, athletes with high physical activity (PA) (>3000 Kcal/day), diabetes, multiple allergies, intestinal or any other dis- eases/conditions that would worsen adherence to the measurements or treatments. The study was conducted at IMIM-Hospital del Mar Medical Research Institute (Barcelona, Spain) from April 2012 to September 2012. Subjects were allocated, by generating random numbers, to one of 3 sequences of administration of 25 mL/day of raw: (a) VOO; (b) func- tional OO enriched with its own PC; and (c) functional OO enriched with its own PC plus complementary phenols from thyme. The statistician was who generated the random allocation sequence while the researcher was who enrolled participants and the doctor who assigned participants to interventions according to the random sequence. Because all participants received each one of the three OOs, restrictions such as blocking were not necessary. The flow-chart of the study was previously described by Pedret et al. [20]. The intervention periods lasted 3 weeks with an ingestion of 25 mL/day raw OO distributed among meals and preceded by 2-week wash-out periods. For the wash-out period, a com- mercial common OO (blend of refined and a small percentage of VOO) kindly provided by Borges Mediterranean Group was used. To avoid an excessive intake of antioxidants and PCs during the clinical trial period, participants were advised to limit the consumption of polyphenol-rich foods (such as vegetables, fruit, coffee, and olives). To register the amount of OO consumed, the volunteers had to return the OO containers to the center after each OO intervention. Participants with less than 80% of treatment adherence (≥5 full OO containers returned) were considered non- compliant subjects for this treatment. PA was evaluated by a questionnaire at baseline and cessation of the study. The trial was conducted in accordance with the Helsinki Declaration and the Good Clinical Practice for Trials on Medical Products in the European Community. All participants provided written informed consent, and the local institutional ethics committees approved the protocol (CEIC-IMAS 2009/3347/I). The protocol was registered with the International Standard Randomized Controlled Trial register (www.controlled- trials.com (accessed on 10 August 2022); ISRCTN77500181) and followed CONSORT- guidelines (Supplementary Table S2). 2.3. Sample Size and Power Analysis A sample size of 30 individuals allows at least 80% power to detect a statistically significant difference among groups of 3 mg/dL of HDL cholesterol, assuming a dropout rate of 15% and a Type I error of 0.05 (2-sided). Antioxidants 2022, 11, 1964 4 of 16 2.4. Anthropometric and Cardiovascular Clinical Measurements Blood pressure (BP), BMI, endothelial function (determined by ischemic reactive hyperaemia (IRH) as previously described [13]), and cardiovascular risk (assessed with the calibrated Framingham function [21]) were measured before and after each intervention. 2.5. Systemic Biomarkers Blood samples were collected following fasting for at least 10 h. Plasma samples were obtained by centrifugation of EDTA whole blood directly after being drawn and were preserved at −80 ◦ C until use. Plasma glucose, total-cholesterol (TC), and triglyceride (TG) levels were measured using standard enzymatic automated methods and apolipoprotein A1 (ApoA1) and apolipoprotein B-100 (APOB100) by immunoturbidimetry in a PENTRA- 400 autoanalyzer (ABX-Horiba Diagnostics, Montpellier, France). HDL-C was measured using an accelerator selective detergent method (ABX-Horiba Diagnostics, Montpellier, France). LDL-C was calculated by the Friedewald equation whenever TGs were inferior to 300 mg/dL. 2.6. Sample Preparation and 1 H-NMR Spectroscopy Serum samples (200 µL) were combined with 400 µL of saline solution (100% D2O), mixed by vortexing, and centrifugated at 13,000× g for 10 min. The supernatant (550 µL) was transferred to a 5 mm internal diameter NMR tube. The metabolic profiles of the serum samples were then characterized by 1 H-NMR spectroscopy using a 500 MHz Bruker NMR spectrometer operating at 310 K. For each serum sample a water-suppressed Carr-Purcell- Meiboom-Gill (CPMG) spin-echo spectrum was acquired using 8 dummy scans followed by 128 scans collected into 64K data points. 2.7. NMR Data Processing 1 H-NMR spectra were manually corrected for phase and baseline distortions. Chemi- cal shifts in the spectra were referenced to the anomeric proton of α-glucose at 5.223 ppm. Spectra were digitized using an in-house Matlab (version R2009b, The Mathworks, Inc.; Natwick, MA, USA) script. Resonances arising from imperfect water saturation were removed to minimize distortions to the spectral baseline. 2.8. NMR Data Analysis 2.8.1. Multilevel-OPLS-DA Analyses Metabolomics data from human studies are characterized by large variations between the subjects. Therefore, subtle treatment effects can be easily overlooked due to the con- siderable variation among participants due to age, disease state, genetics and the like. In addition, the impact of the treatment effect may differ between the subjects. Cross-over designed experiments are particularly useful to tackle these problems, since each subject in the study act as their own controls, making the data paired. A specific limitation of using typical techniques for metabolomics studies such as PCA, PLS or OPLS in a crossover experiments is that the net treatment effect is not separated from the biological variation between the subjects. Therefore, to identify metabolic differences between treatment groups we used a multilevel orthogonal partial least squares discriminant analysis (M-OPLS-DA) modelling [22,23]. It combines the variation splitting property of multilevel simultaneous component analysis (MSCA) with the analysis of the within subjects variation with an OPLS-DA method. MSCA decomposes the data into an offset term, a between-subjects part and a within- subjects part [24,25]. Decomposition of the variation terms is performed in two consecutive centering steps. The first steps is applied on the entire metabolomic data set X and results in a offset term and a mean-centered data block: X = 1 L xTm + Xc = Xm + Xc Antioxidants 2022, 11, 1964 5 of 16 where L is the total number of spectra, Xc contains the mean-centered data, 1L contains ones, and xTm contains the mean values for each column in X. The matrix Xm is defined as the offset term. Then, a second mean-centering is performed per subject i over the K interventions: Xci = 1K xTbi + Xwi where Xci contains the mean-centered data for subject i, 1K contains ones, xTbi contains the mean for each subject i, and Xwi contains the mean-centered data per subject i. Con- catenating the matrices for each individual, we can write the mean-centered data matrix as: Xc = Xb + Xw where Xb contains the between-subjects variation and Xw contains the within-subjects vari- ation. Therefore, the metabolomic data matrix can be written as three terms containing the offset, the between-subjects variation (biological effect), and the within-subjects variation (treatment effect): X = Xm + Xb + Xw Variation splitting was performed using Matlab with the routines available from van Velzen et al. [22]: http://www.bdagroup.nl/content/Downloads/software/software.php (accessed on 1 June 2018). After variation splitting, OPLS-DA is performed on the within-subject data to find systematic differences among the treatment groups. The within-subject variation contains both variation that is equal for all subjects as well as variation different between subjects. An OPLS-DA on the within-subjects matrix is used to focus on the similarity in the treatment effect between subjects. The OPLS-DA was introduced as modification of the PLS-DA method that incorporates orthogonal signal correction (OSC) filters to discriminate between two or more groups using multivariate data [26,27]. PLS is a multivariate regression method based on projections that models the relationship between the metabolic data set X and the response variable Y, but also models the structure of X and Y: X = TP0 + E Y = UC0 + F where T and U are the scores matrices, P and C are the loadings matrices, and E and F are the matrices of residuals of X and Y, respectively. Conversely, OPLS separates the variation in T derived from the PLS in two blocks of variation, predictive (Tp ) and orthogonal (To ) variation: X = Tp P0 p + To P0 o + E Y = UC0 + F The predictive block contains the correlated variation between X and Y, whereas the orthogonal block contains the uncorrelated variation between X and Y. The advantage of OPLS compared to PLS is that a single component (the predictive component) is used as a predictor of the response Y, while the other components describe the other variation that is orthogonal to Y. Therefore, the main advantage of OPLS lies on enhancing model interpretation by forcing all Y-related information into a single component. However, it is worth noting that in terms of predictive power, both PLS and OPLS give similar results [26]. OPLS-DA Loading coefficient plots were generated by back-scaling transformation where covariance is plotted between the Y-response matrix and the signal intensity of the metabolites in the NMR data (X). These plots are coloured based on the correlation coefficient (r2 ) between each metabolite and the Y-response variable, with red indicating strong significance and blue indicating weak significance. The predictive performance (Q2 Y) of the model was calculated using a 7-fold cross-validation approach and model validity was established by permutation testing (1000 permutations). Antioxidants 2022, 11, 1964 6 of 16 2.8.2. Machine Learning Analyses Metabolomics datasets are often high-dimensional with non-linear and complex inter- actions among metabolites. Therefore, we further analysed the metabolomics data using machine learning (ML) models based on trees: random forest (RF). Unlike classical linear model-based statistical methods, RF are fully non-parametric model-free methods that capture complex interaction dependency patterns within predictor features affecting the phenotype. In addition, they provide several variable importance measures (VIMs) that can be used to identify relevant features. RF is an ensemble machine learning method based on “growing” hundreds or thousands of decision trees that uses the average output of all the trees to predict the Y-response [28]. Randomness in the trees is introduced by two elements. First, each tree is built from a random subset sample (bootstrap sample), with replacement, of the original data. Second, at each split in the tree building process, only a random subset of predictor features (mtry) is considered from all candidate predictors. Among these variables, the one providing the best split based on a specific criterion is selected. This has the effect of making a forest with a large number of uncorrelated trees. Then, although each individual tree may perform poorly when predicting the outcome and different trees can give different results, the ensemble of trees have better predictive capability. The advantage of the RF is that the observations not used for the construction of a specific tree (termed out-of-bag (OOB) observations) may be used to estimate the VIM. In particular, we used the mean minimal depth, which is calculated based on the position of the features in the decision tree. Thus, unlike other VIM such as the permutation variable importance, it is only based on the structure of the forest and independent of prediction errors [29]. The RF models and the minimal depth distribution were calculated using the “ranger” and “randomForestExplainer” R packages. However, VIMs do not provide the sign of the association with the response variable. Therefore, to facilitate model interpretation, the contribution and effect of each selected feature (i.e., metabolite) in predicting the treatment group was determined by the exact computation of SHapley Additive exPlanations (SHAP) scores by leveraging the internal structure of RF models [30]. The exact computation of SHAP values guarantees that explanations are always consistent and locally accurate. SHAP values determine the importance of a specific value in a specific feature by comparing the model prediction with and without the feature for each individual. Therefore, the same feature with a specific value may have different SHAP valued for different individuals depending on the interactions with other features of that individual. The SHAP scores were calculated and plotted using the R packages “treeshap” and “SHAPforXGBoost”, respectively. A drawback of VIMs in RF is that they are not directly related to the statistical signifi- cance. Therefore, we further applied an all-relevant ML feature selection strategy based on applying RF iteratively as implemented in the Boruta algorithm [31]. The Boruta approach has been recently proposed as one of the best-performing RF-based variable selection methods for high-dimensional omics datasets [32]. It performs variables selection in four steps: (a) Randomization: by creating a duplicate copy of the original features randomly permutate across the observations (called shadow features); (b) Model building: build a RF with the extended data set (original + shadow features) to compute the normalized permu- tation VIM Z-scores for each feature; (c) Statistical testing: find those relevant features with a VIM higher than the shadow feature with the maximum VIM (MZSF) using a Bonferroni corrected two-tailed binomial test (PBonferroni ). Predictor features with significantly higher, significantly lower, or non-significantly different VIM Z-scores than expected at random compared to the MZSF are deemed important (selected), unimportant (rejected), or tenta- tive, respectively; and (d) Iteration: Unimportant and shadow features are removed and the previous steps are repeated until the status of all features is decided or a predefined number of iterations has been performed. All RF models were calculated using 5000 trees and a number of features (ntree) randomly sampled at each split given by the rounded down square root of the number of features (the mtry recommended for classification problems). Antioxidants 2022, 11, 1964 7 of 16 In addition, in the Boruta algorithm we used 500 iterations and a confidence level cut-off of 0.005 for the Bonferroni adjusted p values. 2.9. Statistical Analyses Normality of continuous variables was assessed by normal probability plots. Non- normally distributed variables were log transformed when necessary. Non-compliant volunteers were excluded from the analysis. To compare means (for normal distributed variables) or medians (for non-normal distributed variables) among groups, ANOVA and Kruskal-Wallis test were performed, respectively. The χ2 and exact F-test, as appropriate, were used to compare proportions. A general linear model for repeated measurements was employed to assess the comparisons of changes in intra- (pre-treatment vs. post-treatment) and inter-interventions (FOOT vs. FOO vs. VOO). A value of p-value < 0.05 was consid- ered significant. Spearman correlations were performed and a p < 0.002 was considered significant due to the assessment of 25 variables (Bonferroni p-value). A PLS model was also used to represent the associations between metabolites and clinical cardiovascular measures. Statistical analyses were performed by SPSS 13.0 software (IBM Corp, New York, NY, USA) and Matlab (version R2009b, The Mathworks, Inc.; Natwick, MA, USA). 3. Results 3.1. General Characteristics, Anthropometric Measurements, and Cardiometabolic Parameters Participants’ baseline characteristics are shown in Table 1. Order of intervention was not found to induce any significant effects in any of the variables measured. No changes in daily energy expenditure in leisure-time PA were observed from the beginning to the end of the study. No changes were observed in blood pressure, BMI, and cardiovascular risk throughout the study. IRH improved after FOOT intervention as previously described [13]. No changes were observed in lipid profile and glucose throughout the study. The three OO interventions were well tolerated by all volunteers who did not reported adverse events. Table 1. Baseline characteristics of the participants. Order 1 (n = 11) Order 2 (n = 11) Order 3 (n = 11) p-Value Anthropometric and Cardiovascular Clinical Measurements Sex: man 5 (45.6%) 7 (63.6%) 7 (63.6%) 0.742 Age (years) 54.91 ± 12.57 55.27 ± 11.88 55.45 ± 7.84 0.856 BMI (Kg/m2 ) 25.63 ± 3.68 26.31 ± 5.25 27.85 ± 4.71 0.529 Physical activity (Kcal/week) 3498.75 (1755.00; 8092.50) 1188.75 (742.50; 1687.50) 3322.50 (861.25; 3663.75) 0.094 Ischemic reactive hyperemia 268.95 ± 344.05 60.374 ± 74.63 177.51 ± 174.01 0.159 (IRH) Cardiovascular risk 4.472 ± 2.424 5.525 ± 2.693 4.970 ± 2.124 0.601 Systemic Lipid Profile and Glycaemia Total cholesterol (mg/dL) 228 ± 43 232 ± 33 219 ± 31 0.680 Triglycerides (mg/dL) 94 (75; 149) 119 (95; 168) 117 (81; 126) 0.517 Glucose (mg/dL) 89 ± 12 93 ± 13 91 ± 11 0.683 HDL-cholesterol (mg/dL) 53 ± 13 53 ± 13 53 ± 20 0.992 LDL-cholesterol (mg/dL) 150 ± 32 152 ± 28 142 ± 26 0.700 Apolipoprotein-A1 (g/L) 1.4 ± 0.2 1.4 ± 0.2 1.5 ± 0.2 0.458 Apolipoprotein-B100 (g/L) 1.2 ± 0.2 1.2 ± 0.2 1.1 ± 0.2 0.529 Values expressed as mean ± S.D. or median (25th to 75th percentile). 3.2. Serum Metabolic Phenotyping A significant MSCA model (Q2 Y = 0.86; p < 0.001) was obtained comparing the serum metabolic phenotypes of individuals after receiving the VOO and FOOT interventions (Figure 1a). Compared to VOO, FOOT ingestion resulted in lower circulating amounts of valine, alanine, glutamine, histidine, and glycine, and decreased levels of citrate, creatine, creatinine, glucose, dimethylsulfone, acetoacetate, and dimethylamine (DMA) in serum. Antioxidants 2022, 11, 1964 8 of 16 The model comparing the serum metabolic profiles after FOO and VOO intakes was also significant (Q2 Y = 0.97, p < 0.001) (Figure 1b). Participants supplemented with FOO had higher amounts of 3-hydroxybutyrate (3-HB), acetate, and acetone in serum. In a similar manner to FOOT, FOO reduced levels of valine, alanine, glutamine, histidine, creatine, creatinine, and DMA in serum compared to VOO. In addition, FOO decreased lactate and isoleucine versus VOO. A significant model was also obtained comparing the serum profiles after FOOT and FOO intakes (Q2 Y = 0.97; p < 0.001) (Figure 1c). Serum collected after FOO intake contained higher amounts of acetone, acetoacetate, and Antioxidants 2022, 11, x FOR PEER REVIEW 9 of3-HB 17 and lower amounts of glutamine and alanine compared to that after FOOT intake. a b c Figure 1.1.MSCA Figure MSCA models modelscomparing comparingthe serum metabolic the serum profiles from metabolic studyfrom profiles participants. MSCA study participants. MSCA coefficient plots comparing the metabolic profiles at (a) comparing the metabolic responses to FOOT coefficient plots comparing the metabolic profiles at (a) comparing the metabolic and VOO intakes, (b) comparing the metabolic responses to FOO and VOO intakes, and (c) com- responses to FOOT and VOO paring the intakes, metabolic(b) comparing responses the metabolic to FOOT responses and FOO intakes. DMA, todimethylamine; FOO and VOODMSO2, intakes, and (c) comparing dime- thylsufone; the metabolicFOO, functionalto responses OO enriched FOOT andwith FOOits intakes. own phenolic DMA,compounds; FOOT, functional dimethylamine; DMSO2,OO dimethylsufone; enriched with its own phenolic compounds plus additional complementary one from thyme; VOO, FOO, functional OO enriched with its own phenolic compounds; FOOT, functional OO enriched with virgin olive oil; 3-OHB, 3-hydroxybutyrate. its own phenolic compounds plus additional complementary one from thyme; VOO, virgin olive oil; 3-OHB, 3-hydroxybutyrate. Antioxidants 2022, 11, 1964 9 of 16 The OPLS-DA results were further validated by ML analyses (Figure 2). Glutamine, creatine, creatinine, DMA, and histidine were consistently decreased after treatment with either FOOT (Figure 2a–c) or FOO (Figure 2d–f) compared to VOO. In addition, FOOT decreased the plasma levels of glycine and DMSO2 compared to VOO (Figure 2c), while FOO also decreased the circulating alanine concentrations but increased the plasma levels of acetone and 3-HB compared to VOO (Figure 2f). Glutamine, creatinine and histidine were the most important metabolites differentiating between the treatments with either FOOT or FOO and VOO. Finally, ML models identified three relevant metabolites discriminative of the treatment with the functional OOs (Figure 2g–i). The circulating levels of isoleucine Antioxidants 2022, were11, xincreased FOR PEER REVIEW in the FOOT compared to the FOO, whereas the levels of acetoacetate and 10 of 17 3-HB were increased in the FOO compared to the FOOT. Figure 2. MachineFigure 2. Machine learning learning results results comparing comparing the the metabolites metabolites from from study participants. study participants. (a–c) Distri- (a–c) Distribu- bution of features’ (metabolites) mean minimal depth, SHAP summary plot, and features selected tion of features’ (metabolites) mean minimal depth, SHAP summary plot, and features selected by the by the Boruta algorithm for the comparison FOOT vs. VOO, (d–f) the comparison FOO vs. VOO, Boruta algorithm forandthe comparison (g–i) FOOTFOOT the comparison vs. VOO, (d–f)The vs. FOO. the comparison minimal depthFOO vs. VOO,plot distribution andshows (g–i) the the features with their distribution of minimal depth and its mean, where the importance of a feature increases with decreasing mean values. In the SHAP summary plot, each dot represents and individual sam- ple. The X-axis represents the SHAP value: the impact of a specific feature (metabolite) on the treat- ment group prediction of a specific individual. Features are sorted in decreasing order based on their overall importance for final prediction (average SHAP values shown in bold). Colours repre- Antioxidants 2022, 11, 1964 10 of 16 comparison FOOT vs. FOO. The minimal depth distribution plot shows the features with their distribution of minimal depth and its mean, where the importance of a feature increases with decreasing mean values. In the SHAP summary plot, each dot represents and individual sample. The X-axis represents the SHAP value: the impact of a specific feature (metabolite) on the treatment group prediction of a specific individual. Features are sorted in decreasing order based on their overall importance for final prediction (average SHAP values shown in bold). Colours represent the values of the metabolites normalized concentrations, ranging from yellow (high concentrations of the specific metabolite) to purple (low concentrations of the specific lipid). The Boruta results are shown Antioxidants 2022, 11, x FOR PEER REVIEW as boxplots 11 of 17 of Variable Importance Measure (VIM) for each selected relevant feature. In the boxplots, the red dot represents the mean and the colour bar above each plot indicates the sign of the association among the feature tions of with the treatment the specific group, metabolite) with(low to purple red concentrations and green indicating a decrease of the specific lipid).or increase The Boruta in re-the FOOT vs. sultsVOO, FOO as are shown vs.boxplots VOO, or of FOOT Variablevs. FOO, respectively. Importance Significant Measure (VIM) for eachfeatures selected were identified relevant fea- using ture. trees, 5000 In the 500 boxplots, the red iterations, dot and represents