Statistical analysis

Biochemical, qPCR and phenotypic data were analysed using and graphpad software. Differential effects were assessed on log2 transformed data by performing ANOVA followed by Sidak post-hoc tests. p-values < 0.05 were considered significant.

Hierarchical clustering of lipid quantification data was performed using R (R Development Core Team, 2018) with the heatmap.2 function from the package, gplots. Data were log2 transformed, then centred and scaled by lipid. Hierarchical clustering was applied to the samples and the lipids using 1-Pearson correlation coefficient as distance and Ward’s criterion (Ward.D2) for agglomeration. All the data represented on the heat map had adjusted p-values <0.05 for one or more comparisons performed with an analysis of variance.

Microarray data were analyzed using R and Bioconductor packages (www.bioconductor.org, v 3.0), as described in GEO accession {"type":"entrez-geo","attrs":{"text":"GSE123354","term_id":"123354"}}GSE123354. Raw data (median signal intensity) were filtered, log2 transformed, corrected for batch effects (microarray washing bath) and normalized using quantile method⁶⁴.

A model was fitted using the limma lmFit function⁶⁵ considering array weights using arrayWeights function. Pair-wise comparisons between biological conditions were applied using specific contrasts. A correction for multiple testing was applied using Benjamini-Hochberg procedure⁶⁶ for False Discovery Rate (FDR). Probes with FDR ≤ 0.05 were considered to be differentially expressed between conditions.

In addition to the differential analysis, a multivariate exploratory analysis was performed. A Sparse Partial Least Squares Discriminant Analysis⁶⁷ (sPLS-DA) was conducted using mixOmics package⁶⁸ under to select the most discriminative variables (genes) that help classify the samples according to their experimental conditions among their expression values. Twenty iterations of 5-fold cross-validation was used to evaluate the model performance for the selection of the most informative components (6 components chosen) using all the variables (PLS-DA). Then a “sparse” PLS-DA model was parametrized selecting the first 100, 120, 80, 40, 100 and 20 (chosen according the performance results of 20 iterations of 5-fold cross-validations) most discriminant variables on the components 1 to 6 respectively.

Hierarchical clustering was applied to the samples and the differentially expressed probes using 1-Pearson correlation coefficient as distance and Ward’s criterion for agglomeration. The clustering results are illustrated as a heatmap of expression signals. Gene network and enrichment of KEGG pathways was either performed using the online software STRING V.11⁶⁹ or Metascape⁷⁰. Correlation graphic chart was generated using the chart correlation function from the Performance Analytics package.