Phospholipidome of endothelial cells shows a different adaptation response upon oxidative, glycative and lipoxidative stress

Endothelial dysfunction has been widely associated with oxidative stress, glucotoxicity and lipotoxicity and underlies the development of cardiovascular diseases (CVDs), atherosclerosis and diabetes. In such pathological conditions, lipids are emerging as mediators of signalling pathways evoking key cellular responses as expression of proinflammatory genes, proliferation and apoptosis. Hence, the assessment of lipid profiles in endothelial cells (EC) can provide valuable information on the molecular alterations underlying CVDs, atherosclerosis and diabetes. We performed a lipidomic approach based on hydrophilic interaction liquid chromatography-tandem mass spectrometry (HILIC-MS/MS) for the analysis of the phospholipidome of bovine aortic EC (BAEC) exposed to oxidative (H2O2), glycative (glucose), or lipoxidative (4-hydroxynonenal, HNE) stress. The phospholipid (PL) profile was evaluated for the classes PC, PE, PS, PG, PI, SM, LPC and CL. H2O2 induced a more acute adaptation of the PL profile than glucose or HNE. Unsaturated PL molecular species were up-regulated after 24 h incubation with H2O2, while an opposite trend was observed in glucose- and HNE-treated cells. This study compared, for the first time, the adaptation of the phospholipidome of BAEC upon different induced biochemical stresses. Although further biological studies will be necessary, our results unveil specific lipid signatures in response to characteristic types of stress.

Our results show for the first time that the lipidome of EC is exquisitely responsive to diverse stress conditions, and thus, may mediate specific adaptive changes.

Results
Lipids have recently emerged as key mediators in the onset of chronic inflammatory diseases characterized by endothelial dysfunction and oxidative stress. Hence, we employed a HILIC-LC-MS/MS platform to analyse the phospholipid profile of BAEC treated in control conditions and in response to different stressing agents (H 2 O 2 , glucose, or HNE). The data sets resulting from the HILIC-LC-MS/MS analysis of four sample groups (control, H 2 O 2 , glucose, and HNE) were later subjected to both univariate and multivariate statistical analysis, aiming to identify significant changes occurring in the BAEC phospholipid profile upon induced biochemical stress.
We performed the identification and the relative quantification of PL species belonging to 8 different classes: phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylserine (PS), phosphatidylglycerol (PG), phosphatidylinositol (PI), lyso-PC (LPC), cardiolipin (CL) and sphingomyelin (SM). The whole list of the 109 PL species (correspondent to the most abundant species in all the identified classes) that were identified and quantified after MS and MS/MS analysis of each sample can be found in Supplementary Table S1. The total chain length (C) and degree of unsaturation (N) are included. Also, the different isomers of the same classes that bear different esterified fatty acids and correspond to each C:N composition were included. These isomers cannot be resolved by the LC-HILIC method, but the fatty acyl composition was determined by MS/MS analysis. Negative ion mode MS/MS data were used to analyse fatty acid carboxylate anions fragments, which allowed to assign the fatty acyl chains esterified to the PL molecular species. For the relative quantification of all the PL listed in Supplementary Table S1, the peak areas of the extracted ion chromatograms (XICs) of each PL species (C:N) within each class were normalized using the peak area of the internal standard (IS) selected for the class. Data were subsequently autoscaled and then subjected to a principal component analysis (PCA) to display the clustering trends of the four experimental groups of BAEC: control, H 2 O 2 -, glucose-, and HNE-treated. The PCA showed that all the groups were separated from each other in a two-dimensional score plot which represented the analyses describing 66.8% of the total variance, including principal component 1 (51.1%) and principal component 2 (15.7%), where principal component 1 was the major discriminating component (Fig. 2).  Additionally, we carried out a hierarchical clustering analysis (HCA) on the phospholipid data sets from the four conditions (Fig. 3). The resulting HCA dendrogram depicted a noticeable separation of the four data sets of control, H 2 O 2 , glucose, and HNE. The first level of separation was evidenced between H 2 O 2 -treated samples and the remaining conditions. The second level of separation distinguished control samples from the two remaining groups (glucose and HNE-treated). However, glucose-and HNE-treated cells differentiated in two clusters in the third level of separation.
Then, we performed a projection to latent structures discriminant analysis (PLS-DA) in order to maximize the phenotypic classification of samples, which showed the performance statistics of R2X = 0.94241, R2Y = 0.96805 and a high prediction parameter Q2 of 0.80026 (X) and 0.90986 (Y). The four groups were well separated in the resulting two-dimensional score plot (Fig. 4). The PLS-DA score plot described 65.8% of the total variance, including component 1 (16.8%) and component 2 (49%). Along with component 2, control samples were scattered at the central region of the plot. Glucose-and HNE-treated samples were scattered on the top region of the plot. H 2 O 2 -treated samples formed the only group that was scattered at the bottom region of the plot. Along with component 1, control-and glucose-treated samples were scattered at the left side of the plot, while HNE-and H 2 O 2 -treated samples were scattered on the right side of the plot.
Besides multivariate statistical analyses, we used another approach to facilitate the interpretation of the extensive dataset produced in the present study, namely the analysis of semi-quantitative phospholipidomic features across distinct depths of detail.
For a clearer overview of the phospholipidome, we included the analysis of PL by features, assessing the adaptation of all the PL species bearing the same number of unsaturation, or bearing the same total carbon chain length, as already performed by other authors 25 . Therefore, the cumulative levels of all the PL species comprised of the same number of unsaturation, ranging from zero, one, two, three and four double bonds, were summarized in the singular lipidomic features PL-DB0, PL-DB1, PL-DB2, PL-DB3, and PL-DB4, respectively. When comparing globally the degree of unsaturation observed for BAEC PL upon the four tested conditions (control, H 2 O 2 , glucose and HNE), the highest cumulative levels of (poly-)unsaturated species (PL-DB1, PL-DB2, PL-DB3, and PL-DB4) were always observed after the treatment with H 2 O 2 (Fig. 5). Conversely, treatment of BAEC with glucose and HNE downregulated the levels of PL-DB1, PL-DB2, PL-DB3 and PL-DB4 when compared with control (Fig. 5).
Only when comparing HNE with control, changes in PL-DB4 were not statistically significant.
Analogously, the cumulative levels of all the PL species comprised of the same number of carbons in their hydrocarbon chain were summarized in the singular lipidomic features PL-C32, PL-C34, PL-C36, PL-C38, and PL-C40, respectively. Regardless of the treatment, PL-C36 species were always the most abundant in BAEC. Treatment of BAEC with H 2 O 2 lead to an increase of the features PL-C34, PL-C36, and PL-C38 when compared with control (Fig. 5). Conversely, we observed a downregulation of PL-C32, PL-C34, PL-C36, PL-C38, and PL-C40 for BAEC treated with glucose, and a downregulation PL-C32, PL-C34, PL-C36, and PL-C40 for BAEC treated with HNE, in comparison with control (Fig. 5). We did not observe any statistically significant variation for PL-C32 and PL-C40 after treatment with H 2 O 2 , nor for PL-C38 after treatment with HNE.
For a more detailed interpretation of the data, we further addressed the adaptation of single PL species induced by the different stressing treatments. The most abundant PC molecular species was PC (34:1) followed by PC (36:2), in all the conditions. We observed a statistically significant increase of the levels of PC (34:1), PC  Table 1).
We did not observe any statistically significant difference for the levels of the molecular species belonging to the CL and LPC classes.

Discussion
Lipids have been highlighted as biomolecules involved in the onset of CVDs, a family of chronic inflammatory diseases in which vascular pathobiology is associated with endothelial dysfunction. Cellular lipid profiling can provide evidence at the molecular level for the role of lipids in these diseases. In our study, we performed a phospholipidomic profiling on BAEC subjected to H 2 O 2 , glucose, or HNE treatment, aiming to elucidate the adaptations in the PL of EC upon three models of biochemical stress associated with the onset and/or progression of CVDs.
In a first instance, the multivariate analyses that we performed (PCA, PLS-DA, and HCA) indicated that the phospholipid profiles of BAEC were significantly altered in response to each one of the induced biochemical stresses (H 2 O 2 , glucose or HNE). The two-dimensional plots of PCA (Fig. 2) and PLS-DA (Fig. 4), along with the dendrogram of HCA (Fig. 3), showed that the most evident differentiation of the BAEC-treated phospholipidome, in comparison with controls, occurred in the model of oxidative stress, 24 h after exposure to 1 mM H 2 O 2 . This suggests that the H 2 O 2 -mediated toxicity induced more particular and specific changes in the PL profile of BAEC when compared with glucotoxicity and lipotoxicity, presumably due to the adaptation mechanisms established by the cells to survive the exposure to oxidative stress. Yang and co-authors 24 studied a hybrid cell line (EA. hy296) subjected to oxidative stress, and using phospholipidomic data found that samples were clearly clustered depending on the different exposure times to 0.2 mM H 2 O 2 (0, 1, 2, 3 and 6 h).
Besides the present study, no other published works have reported the adaptation of the phospholipid profile of EC upon hyperglycemia. In our study, BAEC treated with 25 mM glucose for 24 h showed a significant alteration of the phospholipidome. PCA, HCA and PLS-DA (Figs 2, 3 and 4, respectively) all exhibited good clustering of the phospholipid profile of glucose-treated samples when compared to controls, differing from the trend observed for H 2 O 2 -treated samples. The clustering evidenced by HCA, in which high glucose samples are relatively closer to controls in comparison to H 2 O 2 samples, suggests that the cellular effects mediated by glucose affected the phospholipid profile of BAEC via less harsh changes when compared to a directly induced acute oxidative stress.
In this study, we also assessed the phospholipidome of BAEC upon HNE exposure. Although the biological effects of HNE on several EC types have already been widely studied (as reviewed by Chapple and co-authors 12 ), this is also the first study that aimed to evaluate the variations in the phospholipid profile of EC treated with an ALPP. In this study, we used 10 µM HNE, a concentration in the range of those measured in disease states 12 . BAEC treated for 24 h with HNE showed a significant phenotypical differentiation of their phospholipid profile, as shown in Figs 2, 3 and 4.
Treatment of BAEC with H 2 O 2 increased the cumulative levels of (poly-)unsaturated PL species PL-DB1, PL-DB2, PL-DB3, and PL-DB4, in comparison with control (Fig. 5). Remarkably, it had been previously described that direct oxidative stress (H 2 O 2 1 mM) induced PL peroxidation within the first hours of treatment 26 , while along 24 hours of treatment, BAEC may have established an adaptive mechanism of re-acylation of unsaturated fatty acyl chains in the PL pool, aimed to counteract the ROS-mediated peroxidation. An increase in (poly)-unsaturated PL species was also reported by Peterson and co-authors 27 , in primary neurocortical cells incubated with H 2 O 2 for 24 h. Hence, this is in agreement with our results that show augmented levels of (poly)-unsaturated species in H 2 O 2 -exposed BAEC. The amount of (poly)-unsaturated PL species in the plasma membrane is known to be an effective modulator of its physical properties, as phase transition temperature 28 . In our study, the up-regulation observed for the (poly)-unsaturated species of BAEC subjected to oxidative stress might have resulted in increased membrane fluidity. We found this adaptation of interest since an augmented membrane fluidity was also observed by Sergent 29 in rat primary hepatocytes subjected to ethanol-induced oxidative stress, which would be by our conjectures. However, it is important to remark that in the present study the phospholipidome composition was not strictly studied for the plasma membrane. Hence PL (poly)-unsaturated species might also belong to other cellular organelles. Additionally, other lipids that were not analysed in the present study (e.g., cholesterol) are also regulators of membrane fluidity.
The attempt to evaluate significant changes in both high and low abundant PL molecular species led us to perform a univariate statistical analysis using autoscaled data. Importantly, whereas highly abundant PL species are central for maintaining structural and biophysical properties, lower abundant PL species might more likely act as mediators of signalling functions.
Treatment with H 2 O 2 resulted in a decrease of two PC molecular species esterified to saturated fatty acids, namely PC (30:0) and PC (32:0), along with an increase of three (poly)-unsaturated PC molecular species, namely PC (34:1), PC (34:2) and PC (36:4) ( Fig. 6 and Table 1 32 , hence the increase of singular unsaturated molecular species correlates with the increase that we observed for PL-DB1, PL-DB2, PL-DB3, and PL-DB4. Due to the relatively small head group and the high unsaturation degree characterizing PE molecules, significant changes in the composition of PE species may affect the curvature of the cell membrane and contribute to an increase in its fluidity 33,34 Of interest, Hailey and co-authors 35 found mitochondrial PE to be the recruiters of autophagy markers in starved mammalian cells, while Rockenfeller and co-authors 36 highlighted that increased levels of PE enhanced the lifespan of cultured mammalian cells via promotion of autophagy. Since high levels of H 2 O 2 are known to induce apoptosis in cultured EC 26 , we postulate that the up-regulation of PE species herein observed is an adaptive mechanism of survival adopted by the cells in response to the early apoptotic signals triggered by the H 2 O 2 treatment. Curiously, H 2 O 2 treatment led to a complex regulation of (poly)-unsaturated PS species in BAEC, since PS (38:4) and PS (40:7) were up-regulated, while PS (36:2) and PS (40:5) were down-regulated. The anti-inflammatory properties of PS in mammalian cells have recently been highlighted 37 , hence the increase of PS (38:4) and PS (40:7) could represent a mechanism promoted by BAEC to compensate the inflammatory state triggered by H 2 O 2 26 . Since the formation of ox-PS in mammalian cell membranes is a hallmark of apoptosis 38,39 , we interpret that the decrease in (poly)-unsaturated PS is a consequence of the early apoptotic signals that might have been triggered by incubation with H 2 O 2 .
Five SM species were up-regulated after treatment with H 2 O 2 , namely SM (34:2), SM (36:2), SM (38:2) and SM (40:2). In eukaryotic cells, SM is known to co-localize and associate with sterol lipids and display structural roles related to cholesterol homeostasis and membrane distribution 40 . The interaction between SM and cholesterol forms ordered domains that are primarily involved in the regulation of cell membrane fluidity. In this sense, the oxidative stress-induced up-regulation of SM species observed in this study might reflect a mechanism of perturbation in the fluidity of the plasma membrane of treated BAEC. Nevertheless, there are other SM cellular functions that deserve discussion. Yang et al. 41 treated yeast cells with increasing H 2 O 2 concentrations and found that cell death effects were strongly abrogated in cells expressing sphingomyelin synthase 1 (SMS1), the principal enzyme catalysing the synthesis of cellular SM species. Since BAEC always maintained cell viability along the 24 h treatment with H 2 O 2 (no cell death was observed), the significant increase in SM molecular species might represent an augmented SMS1 activity, adopted by the cells to abolish the cell death signals induced by H 2 O 2 .
Altogether, the significant variations herein observed for several PL and SM molecular species might arise from a wide array of BAEC adaptive mechanisms occurring upon H 2 O 2 -induced stress, which include response to inflammation, impaired membrane fluidity, and cell death. However, it is important to remark that none of the SCieNtifiC RePORtS | (2018) 8:12365 | DOI:10.1038/s41598-018-30695-0 above-listed responses was measured in the present study, hence future works will be necessary to validate what has been herein conjectured.
The global observation of PL features highlighted that the treatment with glucose reduced the levels of PL-DB1, PL-DB2, PL-DB3, and PL-DB4 when compared with control (Fig. 5). A similar trend was observed for PL-C32, PL-C34, PL-C36, and PL-C38 (Fig. 5). To the best of our knowledge, the effects of hyperglycemia on the phospholipidome of EC have never been reported before. In this regard, it is interesting to pinpoint the experiments of Mīinea et al. 42 , which found a diminished Δ5 desaturase activity and a markedly decreased biosynthesis of arachidonic acid-containing PL in Schwann cells grown in 30 mM glucose.
The majority of the PL species that showed significant variation, reported in Fig. 6 and Table 1, were down-regulated after the treatment with glucose. The alterations induced by high glucose may reflect several mechanisms including biological signalling, redox imbalance and formation of advanced glycation end products (AGEs). Hempel and co-authors 43 observed an increased EC permeability mediated by the activation of PKC, upon hyperglycemia. Duffy and co-authors reported that 72 h hyperglycemia induced apoptosis in human aortic EC (HAEC) 44 . Noteworthy, an excess of glucose can lead to NADPH consumption and formation of sorbitol, which in turn can induce oxidative and osmotic stress, respectively 45 . Several authors have reported the promotion of ROS generation and the redox imbalance mediated by high glucose in EC [46][47][48] . The ability of high glucose to induce 12/15-lipoxygenase (12/15-LOX) in mesangial cells 49 and in EC 50 was also highlighted. Glucose is also known to modify proteins 51 and aminophospholipids (PE and PS) 52 . Altogether, these mechanisms are referred as glucotoxicity 53 and contribute to the endothelial dysfunction that occurs in diabetes 54 . However, the cellular mechanisms linking glucotoxicity and PL turnover in mammalian cells are still unclear. On one hand, the formation of aminophospholipid-AGEs could have contributed to the decreased levels observed for PE and PS species; additionally, PE AGEs were found to trigger oxidative stress 55 , which can be a further explanation of the down-regulation of (poly)-unsaturated PL species observed in glucose-treated BAEC. On the other hand, the decrease of (poly)-unsaturated PL species may have been induced by the radical-based oxidative stress promoted by the 24 h exposition to high glucose. However, the ability of 12/15-LOX to oxidize (poly)-unsaturated fatty acid esterified to PL has recently been suggested 56 , hence a glucose-induced up-regulation of 12/15-LOX could be another explanation for the decrease in (poly)-unsaturated PL.
We analysed the global changes occurring in the PL features of BAEC after the treatment with HNE. The levels of PL-DB1, PL-DB2, and PL-DB3, along with the levels of PL-C32, PL-C34, PL-C36, and PL-C40 were downregulated when compared with controls (Fig. 5). Our work constitutes the first report on the effects of HNE on the phospholipidome of BAEC. The treatment with HNE also led to a down-regulation of all the PL molecular species reported in Fig. 6 and Table 1. HNE induces several intracellular cascades in EC, as reviewed by Chapple and co-authors 12 . Increased levels of this ALPP (1-100 µM) were found in disease states, mediating cellular-damaging pathways, including increased ROS generation 57 and endoplasmic reticulum (ER) stress, altogether leading to EC dysfunction. Although the interplay between HNE-mediated toxicity (lipotoxicity) and the alteration of EC lipidome is mostly unknown at present, the ER stress induced by HNE in EC 58 should be carefully considered, since this organelle is directly involved in the biosynthesis of several PL classes 32 . The down-regulation that we observed for all the PL species can be due to a HNE-induced perturbation of the PL biosynthetic pathways in the ER. However, BAEC treated with HNE were also found to suffer increased intracellular ROS levels, increased apoptosis and down-regulation of antioxidant defences 57 . This oxidative stress response triggered by HNE in BAEC can finally lead to membrane PL peroxidation, which would explain the decreased levels of (poly)-unsaturated PL species that were herein observed. Moreover, similarly to glucose, HNE is able to form covalent adducts with nucleophilic sites in aminophospholipids 59,60 , which could further explain the decrease of PS and PE molecular species in HNE-treated BAEC.
In summary, the results from the present study point out that the phospholipidome of BAEC suffers statistically significant changes upon different biochemical stresses (H 2 O 2 , glucose, and HNE). For the first time, the phopsholipidomic profiling of BAEC was compared between homeostasis, oxidative stress, glucotoxicity and lipotoxicity, and specific lipidome alterations were reported for all the tested conditions. The molecular adaptation observed for the treatment with H 2 O 2 was distinctive when compared with the alterations induced by glucose or HNE. More deeply, we found H 2 O 2 to increase the cellular levels of (poly)-unsaturated PL molecular species, while the same species were down-regulated after the treatment with glucose or HNE. These evidences highlight that PL are fundamental players in the response of vascular cells to such external stresses. A specific adaptation of the whole PL profile of EC can represent a cellular hallmark for the onset and the development of CVDs, atherosclerosis, and diabetes. Nevertheless, the biochemical mechanisms adopted by BAEC that lead to the alteration of the phospholipidome are still unclear, particularly in the case of glucolipotoxicity, and a considerable gap of knowledge exists regarding the lipid profiling of EC in non-homeostatic or pathological conditions. These results open new avenues of research towards further studies necessary to unveil the interplay between biochemical stress, PL turnover, and onset of cardiovascular pathologies and diabetic complications. purification. All the other reagents and chemicals used were of the highest grade of purity commercially available. The water was of Milli-Q purity (Synergy1, Millipore Corporation, Billerica, MA).

Reagents
Cell culture and treatments. Bovine aortic endothelial cells (BAEC) were obtained from Lonza, Inc., (Walkersville, MD) and cultured in RPMI1640 medium supplemented with antibiotics (100 U/mL penicillin and 100 µg/mL streptomycin) and 10% calf serum from Gibco (Life Technologies). For each experiment, BAEC were used between passages 8 and 16 and were grown to near confluence for experiments (80-90% density). Milli-Q water was used as a solvent throughout the experiments (H 2 O 2 , glucose, HNE). Treatment of BAEC with the different stressing agents (H 2 O 2 , glucose, HNE) was carried out in serum-free medium. This induces a near-quiescent state and does not reduce cell viability 61 . Control cells received an equivalent amount of water as required. Cells were subjected to biochemical stress with 1 mM H 2 O 2 , 25 mM glucose, or 10 µM HNE during 24 h. After treatment, cells were collected by scraping in PBS on ice and centrifuged at 1000 rpm for 5 min. Cell pellets were stored at −80 °C. The whole experimental procedure, including cell culture and treatments, lipid extraction and analysis, was repeated four times.
Lipid extraction and quantification of phospholipid content. Thereafter, total lipids were extracted using the Bligh and Dyer method 62 , and the phospholipid (PL) amount in each lipid extract was quantified with the phosphorus assay, performed according to Bartlett and Lewis 63 . For the detailed experimental procedures of lipid extraction and PL quantification, the reader is referred to a previously published work in which the same methodologies were applied 64 . Data and statistical analysis. Phospholipid peak integration and assignments were performed using MZmine 2.30 65 . The software allows for filtering and smoothing, peak detection, peak processing, and assignment against an in-house database. During the processing of the raw data acquired in full MS mode, all the peaks with raw intensity lower than 1e4 were excluded. For all assignments, peaks within 6 ppm of the lipid exact mass were considered. Consequently, assigned PL were further validated by manual analysis of the MS/MS data (Supplementary Table S1). Analysis of the MS/MS spectra acquired in positive ion mode was performed to confirm the identity of the molecular species belonging to the PC, LPC, SM, and PE classes. The fragment ion at m/z 184.07, when observed in the MS/MS acquired in the positive mode, indicates the presence of phosphocholine head group, allowing to identify the precursor ions as belonging to the PC, LPC and SM classes. The species belonging to these three classes were further differentiated by the characteristic retention times. The neutral loss of 141 Da, when observed in the MS/MS acquired in the positive mode, allows identifying precursor ions of the PE class. The MS/MS spectra acquired in negative ion mode were used to confirm the identity of PG, PI and PS phospholipids. In particular, the fragment ion at m/z 171.01 was used to confirm [M-H] − PG class precursor ions and the fragment ion at m/z 241.01 was used to confirm [M − H] − PI class precursor ions. The neutral loss of 87 Da, when observed in the MS/MS spectra acquired in the negative mode, allows identifying precursor ions belonging to the PS class. Negative ion mode MS/MS data were used to identify the fatty acid carboxylate anions fragments RCOO − , which allowed the assignment of the fatty acyl chains esterified to the PL precursor. All the MS/MS fragmentation patterns characteristic for the lipid classes analysed in the present study, acquired in positive and negative ion modes, are available online as supplementary information.

HPLC-ESI-MS and MS
Relative quantitation was performed by exporting integrated peak areas values into a computer spreadsheet (Excel, Microsoft, Redmond, WA). For normalization of the data, the peak areas of the XICs of the PL precursors of each class (listed in Supplementary Table S1), (C:N composition), were divided for the peak area of the IS selected for the class.
SCieNtifiC RePORtS | (2018) 8:12365 | DOI:10.1038/s41598-018-30695-0 Initial data assessment was performed using Metaboanalyst 66 . Shapiro-Wilk test was used for analysing the normality of the BAEC treatment data. Univariate statistical analysis was performed using ANOVA test following post hoc least significant difference (LSD) test, with Benjamin-Hochberg correction for the false discovery rate (FDR). Univariate data processing and analyses were performed using the SPSS software package (IBM SPSS Statistics Version 24). Principal component analysis (PCA), projection to latent structures discriminant analysis (PLS-DA) and hierarchical clustering were performed on auto-scaled data using the R software (R version 3.4.2 downloaded from https://www.R-project.org) with the packages RFmarkerDetector, FactoMineR 67 , Factoextra 68 and Ropls 69 . PCA was performed using the FactoMineR and ellipses were drawn assuming a multivariate normal distribution. Hierarchical clustering was performed using Ward's method using the method hclust. Boxplots were created using the R package ggplot 70 . Bar graphs were created using the software GraphPad Prism 7.
Data availability. The datasets generated during and/or analysed during the current study are available online as supplementary material.