## 2 Materials and Methods
### Multivariate Data Analysis (MVDA)  
The acquired high-resolution GC-TOF-MS datasets were converted to NetCDF files using LECO ChromaTOF-HRT software, version 5.0. These were exported to the XCMS online statistical package ([[https://xcmsonline.scripps.edu](https://xcmsonline.scripps.edu)]([https://xcmsonline.scripps.edu](https://xcmsonline.scripps.edu)), accessed on 10 November 2019), an automated, web-based metabolomics data processing software for peak picking and peak alignment. The parameters for the method were selected for GC-TOF-MS specificities and were as follows: (i) feature detection was set as the centWave method, minimum peakwidth = 0.6, maximum peak width = 1.0; (ii) Rt correction was set as the Obiwarp method, Profstep = 1, alignment was set as _m/z_ width = 0.015, minfraction = 0.5 and bw = 5. The resulting peak list showed 7550 variables (in a combined list) with corrected peak retention times (Rts, min), mass-to-charge ratios (_m/z_) and integrated peak areas.  
The resultant data matrixes obtained from the XCMS statistical package were imported into SIMCA-P software, version 15.0 ([PERSON], Umea, Sweden) and Pareto-scaling was applied for the multivariate statistical analysis. Firstly, to reduce the dimensionality of the data and to summarize the information contained in the datasets, a principal component analysis (PCA, an unsupervised method for MVDA) was used. Here, the [PERSON]'s T-squared (T\({}^{2}\)) test was applied to test for significant differences between the mean vectors (multivariate means) of the multivariate data sets. Secondly, to complement the PCA, a hierarchical cluster analysis (HCA) was used to assess the trends and patterns observed on the PCA [39; 40; 41].  
Orthogonal partial least squares discriminant analysis (OPLS-DA), a supervised MVDA method, was used to extract the maximum amount of information regarding the significant variables from the datasets. This procedure allows for the removal of systematic variation from the experimental data (X variables) that is not correlated to discriminant classes (Y) and [PERSON]'s region, represented by the ellipse and defined by the 95% interval of confidence. The quality of the generated models is described by diagnostic tools for metabolomics, which included: (i) cumulative model variation in the matrix X, (ii) the goodness-of-fit parameter [R\({}^{2}\)X(cum)], (iii) the model variance proportion [Y\({}^{2}\)X(cum)] and (iv) total variation of the matrix X predicted by an extracted component [Q\({}^{2}\)(cum)] [42].  
In addition, the OPLS-DA models were statistically validated using a cross-validated analysis of variance (CV-ANOVA), with a \(p\leq 0.05\) indicating a good model [43], using the SIMCA inbuilt 7-fold (default) CV method and a permutation test (50 permutations). OPLS-DA S-plots were constructed for subsequent interpretation. Based on the S-plots, significant _m/z_ ions with a correlation of [(_p(corr)_) \(\geq\) 10.51 and a covariance of (_p1_) \(\geq\) 10.051 were selected. These ions were further validated using 'variable importance in projection' (VIP) scores of >1 before selection for annotation. Statistically significant features were annotated (tentatively identified) to level 2, as defined by the metabolomics standards initiative [44], based on their mass spectral information using the National Institute of Standards and Technology (NIST) Mainlib [45], and Fiehn metabolomics and Mass Spectrum libraries [46], selecting compounds with a similarity of 75% and above.