**Lipid Annotation by Combination of UHPLC-HRMS (MS), Molecular Networking, and Retention Time Prediction: Application to a Lipidomic Study of in Vitro Models of Dry Eye Disease**

**Romain Magny 1,2, Anne Regazzetti 2, Karima Kessal 1,3, Gregory Genta-Jouve 2,4, Christophe Baudouin 1,3,5, Stéphane Mélik-Parsadaniantz 1, Françoise Brignole-Baudouin 1,2,3, Olivier Laprévote 2,6 and Nicolas Auzeil 2,\***


Received: 10 April 2020; Accepted: 25 May 2020; Published: 29 May 2020

**Abstract:** Annotation of lipids in untargeted lipidomic analysis remains challenging and a systematic approach needs to be developed to organize important datasets with the help of bioinformatic tools. For this purpose, we combined tandem mass spectrometry-based molecular networking with retention time (tR) prediction to annotate phospholipid and sphingolipid species. Sixty-five standard compounds were used to establish the fragmentation rules of each lipid class studied and to define the parameters governing their chromatographic behavior. Molecular networks (MNs) were generated through the GNPS platform using a lipid standards mixture and applied to lipidomic study of an *in vitro* model of dry eye disease, *i.e.*, human corneal epithelial (HCE) cells exposed to hyperosmolarity (HO). These MNs led to the annotation of more than 150 unique phospholipid and sphingolipid species in the HCE cells. This annotation was reinforced by comparing theoretical to experimental tR values. This lipidomic study highlighted changes in 54 lipids following HO exposure of corneal cells, some of them being involved in inflammatory responses. The MN approach coupled to tR prediction thus appears as a suitable and robust tool for the discovery of lipids involved in relevant biological processes.

**Keywords:** lipidomic; liquid chromatography; tandem mass spectrometry; molecular network; dry eye disease; hyperosmolarity

#### **1. Introduction**

Over the last decade, lipids have become a major research topic and are now recognized as key biological compounds displaying various roles in cell functions. They include the coordination of bio-membrane structures, intra- and extra-cellular communication, metabolic efficiency, and signaling

cascades, all of which are critical for cell functionality [1]. Disruption of lipid homeostasis is now recognized to be involved in numerous pathologies such as cancer, diabetes, neurodegenerative disorders or chronic inflammatory diseases [2,3]. Lipids, especially phospholipid and sphingolipid classes, are central in both inflammatory and cell death processes [4,5]. Arachidonic acid, mainly originating from the cleavage of phospholipids, is, for example, widely recognized as a pro-inflammatory fatty acid [6]. Besides, ceramides, a sphingolipid subclass, mediate apoptosis through a caspase-3 dependent mechanism and inflammation through the release of cytokines such as IL-1ß or IL-6 [7].

Nevertheless, lipids encompass a tremendous number of molecular species exhibiting a wide variety of structures. Indeed, the cellular lipidome includes numerous subclasses of sphingolipids, phospholipids, glycerolipids, sterol lipids, and lipid metabolites. Lipidomics, the comprehensive analysis of lipids in biological systems, remains challenging and must involve not only efficient analytical techniques but also appropriate sample processing and integrative computational approaches [8]. Electrospray ionization (ESI) mass spectrometry (MS), either through a shotgun approach or hyphenated to liquid chromatography (LC), has become the gold standard of lipidome study [8,9]. Indeed, lipidomic analysis using an infusion approach represents a useful strategy to easily access a large part of the lipidome [10,11]. Nevertheless, despite fruitful applications, this approach still suffers from ion suppression which limits the analysis of low abundant lipid species [10,12]. In contrast, lipidomic analysis using liquid chromatography hyphenated to mass spectrometry includes a separation step reducing ion suppression and, therefore, improves detection of low abundant lipid species through an increase in sensitivity [9,11,13]. Recent advances highlight that ion mobility, combined to chromatography and mass spectrometry, represents an additional source of information in the case of lipid identification [14,15]. Nevertheless, whatever the benefit of these techniques, the fact is that, in the context of a lipidomic analysis, the large amount of data to be processed requires the development of new data processing approaches, in particular, to simplify the annotation of lipid species while maintaining a high degree of reliability.

In the course of a lipidomics study, the reliable identification of the numerous lipid species detected represents a rigorous and demanding task. Lipids identification is performed taking into account four analytical features: retention time (tR), accurate precursor ion *m*/*z* value, isotopic ratio, and MS/MS data through comparison to reference compounds [16]. In the case of phospholipids, identification deals with the determination of the polar head group and the length of the acyl chains and of their sn1/sn2 location on the glycerol moiety. It must be emphasized that such identification is biologically strongly relevant inasmuch as the nature and position of acyl chains depend on the homeostatic balance between biosynthesis rates, remodeling, and degradation and also reflect the extent of the pool of fatty acids available.

To elucidate the structure of unknown compounds, bioinformatics tools are strongly valuable. Among them, molecular networks (MNs) have recently been proposed, historically in the field of plant secondary metabolites, to identify compounds of biotechnological interest or exhibiting promising pharmacological activity [17–20]. Molecular networks are a computational strategy aimed at organizing and visualizing hundreds of molecules using their MS/MS spectra in accordance to their similarities through the assumption that structurally related molecules display similar product ion spectra [21]. The structural similarity of a set of compounds is visualized in a network which may be generated through online platforms such as GNPS or MetGem [21,22].

Taking into account the LC retention time represents an important benefit for the identification of compounds, as it makes it possible to discriminate isobaric compounds [23]. For this purpose, a pre-processing step of the LC-MS/MS data is required; it may be performed using software such as MzMine 2 [24]. Furthermore, based on structural properties, retention time can easily be predicted with good reliability, and tR values have previously been used to support the identification of numerous metabolite classes, especially lipids [25,26].

The aim of our study was to propose a new approach for the rapid and reliable structural annotation of phospholipids and sphingolipids species using molecular networking and tR prediction. Using 65 commercial lipid standards, we first determined the collision energy conditions required to achieve the fragmentation patterns appropriate to the structural annotation of lipids. To support lipids annotation, commercial standards were further used to define the relationship between lipid structure and tR. Unknown lipids were then identified based on their exact mass measurements and MS/MS fragmentation through molecular networking and tR values.

This approach was used to perform a lipidomic analysis of human corneal epithelial cells exposed to hyperosmolarity (HO)—an *in vitro* model of dry eye disease (DED) [27,28]. Dry eye disease, a chronic multifactorial inflammatory pathology, is characterized by alteration of tear film, cell damage, and inflammation of the ocular surface [29,30]. This very common ocular pathology is also characterized by HO, one of the core mechanisms of DED [31]. Disruption of lipid homeostasis, known to be involved in inflammation and the cell death process, may also be a key feature in the pathophysiology of DED [32]. Thanks to MNs and tR prediction, our lipidomic approach allowed annotation of 150 unique lipid species and highlights homeostasis disruption of 54 lipid species. Several of them are involved in inflammation and cell death.

#### **2. Results and Discussion**

Reliable annotation of unknown compounds using MN highly depends on the quality of the acquired MS/MS spectra [33]. For this purpose, we performed a set of MS/MS experiments for which collision energy was increased step by step in order to optimize the diagnostic fragment ion intensities. Phospholipid annotation needs the presence on the MS/MS spectra of product ions corresponding to the polar head group, the fatty acyl side chains and of the precursor ion. For sphingolipid annotation, fragments corresponding to the sphinganine base moiety and fatty acyl side chain must be detected on MS/MS spectra. Figures 1 and 2 exhibit the main MS characteristics related to PC (16:0/18:1) and Cer (d18:1/16:0), respectively.

#### *2.1. Fragmentation Patterns of Phospholipids*

In the negative ion mode, phosphatidylcholine (PC) are mainly detected as [M−CH3] − ions corresponding to an in-source loss of a methenium and, to a less extent, as formiate ([M+HCOO]−) and acetate ([M+CH3COO]<sup>−</sup>) adducts (Figure 1A). In our study, annotation of a PC was based on the detection of six diagnostic product ion peaks in the MS/MS spectrum of [M−CH3] −.

In the example shown in Figure 1 related to PC (16:0/18:1), the precursor ion was observed at *m*/*z* 744.5540. At low mass, a peak at *m*/*z* 168.0423 corresponded to the deprotonated demethylated phosphocholine ion formed at a 25 eV collision energy (Figure 1B). At higher mass, the peaks at *m*/*z* 255.2334 and *m*/*z* 281.2479 were assigned to oleate and palmitate and exhibited an increased intensity from 20 to 50 eV collision energy (Figure 1C). Two other key product ions at *m*/*z* 480.3098 and *m*/*z* 506.3256 corresponded to demethylated lysophosphatidylcholine LPC (16:0) and LPC (18:1) ions, respectively. They were detected from 20 to 40 eV collision energy with a maximum intensity at 30 eV (Figure 1D). A collision energy ramp between 20 and 40 eV thus appeared to be suitable to obtain the six diagnostic ions with sufficient sensitivity and mass accuracy ( < 10 ppm) (Figure 1E). The ions used to identify 10 standard PC species are compiled in Table S3.

Interestingly, while in the positive ion mode, the phosphocholine product ion at *m*/*z* 184.0733 allows for the highly sensitive detection of PC, and the negative ion mode is essential to perform fatty acyl chains identification [34]. This ionization mode was also successfully applied to identify the fatty acyl chains of other phospholipid subclasses, namely, phosphatidylethanolamine (PE), phosphatidylinositol (PI), phosphatidylserine (PS), phosphatidylglycerol (PG), and phosphatidic acid (PA). It also proved useful to locate the fatty acyl chains on the sn1 and sn2 positions of the glycerol core. Indeed, for PC, PE, and PG species, the intensity of the carboxylate ion peak corresponding to the fatty acid at the sn2 position was always significantly higher than that at sn1. This is in agreement with previously published data [35]. For example, Figure 1C shows an oleate peak at sn2 more intense than the sn1 palmitate in the whole collision energy range. Similarly, the fragment corresponding to demethylated LPC (16:0) formed by the loss of the sn2 oleate from the precursor ion was more intense than the demethylated LPC (18:1) arising from the loss of the sn1 palmitate (Figure 1D). For the lipids belonging to the PS, PI, and PA subclasses, the acyl group at sn1 always led to the more intense product ion peak (Figures S2 and S3) [35]. However, it is noteworthy that the relative intensity of carboxylate product ions displayed in the MS/MS spectra only provided information on fatty acids sn1 and sn2 locations regarding the major regio-isomer. Indeed, the presence in the mixture of a minor amount of the other regio-isomer cannot be excluded. Ensuring it, would need to build a calibration curve using the two pure regio-isomers [36]. In this study, we thus report what is likely to be the major regio-isomer.

**Figure 1.** MS characteristics of PC (16:0/18:1). (**A**) Relative intensities of [M−CH3] −, [M+OAc]−, and [M+HCOO]<sup>−</sup> ions. (**B**–**D**) Fragmentation patterns of the [M−CH3] − ion formed from PC (16:0/18:1) under negative ionization conditions: relative intensities versus collision energy of (**B**) polar head group, (**C**) fatty acyl side chains, and (**D**) demethylated lysophosphatidylcholine ions. (**E**) Tandem mass spectrum of PC (16:0/18:1) [M−CH3] − ion obtained by collision energy ramping from 20 to 40 eV. Blue and red colors relate to the sn1 and sn2 positions, respectively.

The polar head groups were identified owing to the presence of specific fragment ions such as *m*/*z* 168.0428 and *m*/*z* 224.0694 for PC, *m*/*z* 140.0113 and *m*/*z* 196.0380 for PE, *m*/*z* 227.0326 for PG, and *m*/*z* 241.0119 for PI. For PS, an abundant and specific serine loss (87.0326 Da) was observed in the MS/MS spectra of the deprotonated molecules [M−H]<sup>−</sup> together with a glycerophosphate ion at *m*/*z* 152.9958, whereas PA only led to a glycerophosphate ion. Table S3 in the Supplementary Materials lists the diagnostic ions for the 65 phospholipid standard species included in our study.

#### *2.2. Fragmentation Patterns of Sphingolipids*

Ceramides (Cer) are based on a sphinganine backbone amidated by a fatty acyl side chain. Ceramides are also building blocks of more complex sphingolipids such as hexosyl ceramide (HexCer) or sphingomyelins (SM). Under negative ionization conditions, Cer were mainly detected as [M−H]<sup>−</sup> and, to a lesser extent, as acetate ([M+CH3COO]<sup>−</sup>) and formate ([M+HCOO]−) adducts (Figure 2A). As for phospholipids, annotation of individual ceramides was based on the detection of six diagnostic ions which made it possible to determine the double bound number of the sphinganine backbone as well as the FA side chain length. Fatty acyl chain identification was performed thanks to the MS/MS spectra of the [M−H]<sup>−</sup> ion at collision energies of 20 eV and more (Figure 2B,C).

In the example shown at Figure 2, product ions labelled T (*m*/*z* 280.2646), U (*m*/*z* 254.2486) and S (*m*/*z* 296.259) are indicative of the fatty acyl chain (Figure 2B,C). The sphingosine moiety is characterized by the product ions Q at *m*/*z* 263.2379 and P at *m*/*z* 237.2225 formed at collision energies higher than 20 eV (Figure 2D). The [M−H]<sup>−</sup> precursor ion and its product ions [M−H−H2O]<sup>−</sup> and [M−H−2H2O]<sup>−</sup> are detected in the high mass region (Figure 2E). As exemplified by Cer (d18:1/16:0), a collision energy ramping from 20 to 40 eV proved suitable to detect the six diagnostic ions with good sensitivity and mass accuracy ( < 2 ppm), allowing an easy annotation of Cer species. MS/MS spectra of Cer displayed more intense product ion peaks than PC therefore resulting in a better accuracy of *m*/*z* values.

**Figure 2.** MS characteristics of Cer (d18:1/16:0). (**A**) Relative intensities of [M−H]−, [M+OAc]−, and [M+HCOO]<sup>−</sup> ions. (**B**–**D**) Fragmentation patterns of the [M−H]<sup>−</sup> ion formed from Cer (d18:1/16:0) under negative ionization conditions: relative intensities versus collision energy of (**B**) precursor ion, (**C**) fatty acyl side chains, and (**D**) sphinganine ions. (**E**) Tandem mass spectrum of Cer (d18:1/16:0) [M−H]<sup>−</sup> ion obtained by collision energy ramping from 20 to 40 eV. Product ions: T (*m*/*z* 280.2646), U (*m*/*z* 254.2486), and S (*m*/*z* 296.259) are indicative of the fatty acyl chains; Q (*m*/*z* 263.2379) and P (*m*/*z* 237.2225) are indicative of the sphingosine moiety. Product ions are labelled according to Reference [37]. Green and red colors relate to sphingosine and fatty acyl position, respectively.

#### *2.3. Retention Time Prediction*

A typical UHPLC-ESI-MS negative ion mode chromatogram of a mixture of 65 lipid standards representative of the nine studied subclasses is displayed in Figure 3A. In reversed-phase liquid chromatography, an elution of lipids is closely related to the fatty acyl chain lengths, and this property has been widely used in the frame of lipidomic analyses [38,39]. Under our conditions, the FA and lysophospholipids were firstly eluted for 6 min followed by Cer and phospholipids (*i.e.*, PE, PI, PG, and PS) between 6 and 9 min (Figure 3A). Furthermore, in the case of phospholipids, the chromatographic behavior was also dependent on the polar head group, the elution order being for a given fatty acyl chain pattern as follows: PI, PG, PS, PC, and PE (Figure S1). Retention time may thus be considered as a valuable analytical feature helpful in confirming annotation or in highlighting misannotation. However, this requires a robust chromatographic system able to deliver stable and reliable retention times. In our hands, RP-UHPLC operating with a reduced particle size (1.7 μm), associated to a column temperature of 50 ◦C, and optimized elution conditions provided chromatograms with peak widths lower than 20 s and highly reproducible retention times (Table S2).

**Figure 3.** Retention time prediction model based on equivalent carbon number (ECN) of nine subclasses of phospholipids and sphingolipids. (**A**) UHPLC-HRMS chromatogram of a commercial standard lipid mixture in the negative ion mode. (**B**) Theoretical *t*<sup>R</sup> plotted against experimental *t*<sup>R</sup> for standard lipid mixture before and (**C**) after k value fitting. A linear (phospholipids) and polynomial (sphingolipids) regression model was used. Dotted lines represent the 95% confidence interval displaying a relative error which reached 15% on the linear regression graph (B) (k = 2) and did not exceed 5% on the linear regression graph (C). Table: Parameters of relationship between ECN and experimental tR determined using a standard lipid mixture. The general expression of ECN is ECN = NC - k × DB where NC and DB are the total number of carbons and the number of double bonds, respectively. PC: phosphatidylcholine, PE: phosphatidylethanolamine, PG: phosphatidylglycerol, PS: phosphatidylserine, PI: phosphatidylinositol, PA: phosphatidic acid, Cer: ceramide, SM: sphingomyelin, HexCer: hexosyl ceramide.

In order to reinforce lipid annotation using tR, we first used the commercial lipid standard to build a tR predictive model based on the equivalent carbon numbers (ECNs). The general expression of ECNs is ECN = NC – k × DB, where NC and DB are the total number of carbons and the number of double bonds, respectively. For lipid species, k = 2 is usually applied [40].

In the case of phospholipids, especially when containing polyunsaturated fatty acid (PUFA), the linear tR prediction model using k = 2 was fairly poor (Figure 3B). To improve the tR prediction model, we plotted ECN values with experimental tR allowing, next to an appropriate fitting step, to determine more accurate k values. New k values calculated for the different investigated phospholipid subclasses thus improved the linear correlation between ECN and tR (Figure 3C). For example, the difference between the theoretical and experimental tR values for PE (17:0/20:4) was 13% for k = 2 and decreased to 0% using k = 1.25.

In the case of sphingolipids, the linear tR prediction model with k = 2 was suitable for species whose fatty acyl chains did not exceed 22 carbon atoms. However, regarding this lipid class, a polynomial curve led to a better tR prediction model than a linear one. The selection of the polynomial degree was based on the value of the correlation coefficient using a tolerance of 1/100 (*R*<sup>2</sup> > 0.99). For the three investigated sphingolipid subclasses, a quadratic function was finally retained.

Based on this improved fitting step, we predicted more accurately the theoretical tR for all the commercial standards lipids including six phospholipid subclasses and three sphingolipid subclasses (Figure 3C,D). Indeed, the difference between the experimental and theoretical tR values did not exceed 5%, whatever the commercial standard lipid. Consequently, an uncertainty of 5% was retained when tR was used in the purpose to support the annotation of unknown lipids from HCE cell extract and also to discriminate two isobaric lipid species as exemplified later in the text (see Section 2.6).

#### *2.4. Instrument Stability*

The stability of the UHPLC-ESI-MS/MS system was assessed for both MS and chromatography. For this purpose, exact mass measurements were performed on the standard lipid species. High-resolution mass measurements led to a mass accuracy better than 5 ppm for whichever standard lipid species under consideration (Table S2) and the chromatographic system delivered tR values with a deviation within a 3-day period not exceeding 3% (Table S2).

#### *2.5. Lipidic Networking of Human Corneal Epithelial Cells*

Although untargeted analysis by MS constitutes a relevant and powerful tool to characterize a cell lipidome, the annotation of the lipids of interest remains a real challenge. Fortunately, phospholipids and sphingolipids display structural characteristics which make their identification suitable through molecular networking.

In our case, MNs were used as part of a study aimed at assessing the impact on the HCE cell lipidome exposed to HO and were implemented as described hereafter. Tandem mass spectrometry in the DDA mode (see Section 3.3) was used to acquired MS and MS/MS data for 65 commercial standard lipids and for the whole lipids contained in HCE cell extracts. Next to a preprocessing step performed with MzMine 2 (see Section 3.4), preprocessed data were subsequently used to build MNs through the GNPS platform. The MNs thus included three types of nodes corresponding first, to commercial standard lipids, second, to lipids available both as commercial references and identified in HCE cell extracts, and, finally, to lipids only detected in HCE cell extracts (Figures 4 and 5, Table S3).

In such a MN, commercial lipid standards of known structure were thus used to anchor the molecular network within which lipids from HCE cell extracts were clustered according to the similarities of structures that they shared with reference lipids. Thanks to the standard lipids, the key parameters were optimized to provide a reliable and relevant MN. The nodes of the network, corresponding to the MS/MS spectra, were only linked to others if they displayed a common fragmentation pattern, *i.e.*, a minimum number of six identical product ions and/or neutral losses. Moreover, the similarity score between a pair of MS/MS spectra, also called "cosine score" (cos) had to be greater than 0.6. Values selected for the aforementioned parameters were widely used for molecular networking [21].

Figure 4 corresponds to the network of phospholipids; it includes PC, PE, PA, and PG. It was mainly built from the fatty acyl chains located at the sn1 and sn2 positions, and, to a less extent, to the polar head group. For instance, various phospholipid species (*i.e.*, PC, PE, PC-P, PE-P, PG, PA) containing palmitate in the sn1 and sn2 positions were clustered (cluster MN-C1 in Figure 4). The PC (16:0/18:1) and PC (18:0/18:1) were also clustered, as they both contained oleate in the sn2 position and a phosphocholine polar head group (Figure 4). Similarly, PC (18:0/18:1) and PC (18:0/18:2) were clustered, as they contained stearate in the sn1 position (Figure 4). In some cases, depending on the precursor ion intensity, the MN could also include several adducts corresponding to only one PC. This is especially the case for the lipid standard PC (17:0/17:0) observed as [M−CH3] <sup>−</sup>, [M+HCOO]<sup>−</sup> and [M+CH3COO]<sup>−</sup> ion species (cluster MN–C2 in Figure 4). Phospholipids containing ether (PC–O and PE–O) or vinyl ether (PC–P and PE–P) bonds in the sn1 position are displayed on the same MN as the diacylphospholipids. For instance, PE–P (16:0/16:1) included in MN–C1 was connected to PE–P (16:0/16:0), as they both contained a phosphoethanolamine head group and a sn1 palmitoyl moiety. The PE–P (16:0/16:0) was also connected to PC–O (16:0/16:0) and PE (16:0/16:0), as they all contained the same two palmitoyl moieties. In contrast, PE–P (16:0/16:1) was not connected to PE (16:0/16:1), indeed, their MS/MS spectra displayed only four common product ions, this being insufficient to connect them in the MN. To improve data visualization, the MN was organized under two orthogonal axes; the abscissa and the ordinate corresponding, respectively, to the sn1 and sn2 fatty acyl chains of glycerophospholipids. For instance, lipids containing an arachidonate side chain in the sn2 position were displayed on the same line (cluster MN–C3 in Figure 4).

**Figure 4.** Lipidic molecular network of phospholipids including commercial standards and epithelial corneal cell lipids. The MN between (**A**) PC, PE, PA, and PG; (**B**) PS; (**C**) PI; and (**D**) LP. Networking was based on both fatty acyl side chains and polar head groups. The MN was organized along two orthogonal axes: abscissa and ordinate correspond to the sn1 and sn2 fatty acids, respectively. See the text for the explanation of the MN–C1, MN–C2, and MN–C2.

**Figure 5.** Lipidic molecular network of sphingolipids including commercial standards and epithelial corneal cell lipids. MN of (**A**) ceramide, (**B**) sphingomyelin, and (**C**) hexosylceramide. The MN was organized along two orthogonal axes: abscissa and ordinate corresponding to insaturation of the sphinganine base and fatty acid side chain length, respectively.

It must be emphasized that, in contrast to other phospholipids, PS and PI clustering was mainly based on their polar head group because of abundant characteristic fragment ions (Figure 4B,C). The [M−H]<sup>−</sup> precursor ions of PS were dissociated by the loss of the serine moiety leading to an intense [M−H−87.0326]<sup>−</sup> product ion [35] (Figure S2). The MS/MS spectra of the deprotonated PI displayed an inositol phosphate fragment at *m*/*z* 259.0225 and two product ions corresponding to the consecutive loss of one (*m*/*z* 241.0119) and two (*m*/*z* 223.0013) water molecules (Figure S3).

The networks connecting sphingolipids were all organized by subclasses (Figure 5). Ceramides (Figure 5A) were thus clustered separately from sphingomyelins (Figure 5B) and hexosylceramides (Figure 5C). Furthermore, networks of each sphingolipid subclasses displayed clusters depending on sphinganine (d18:0), sphingosine (sphing-4-enine, d18:1) or a sphingadienine (sphing-4,14-dienine, d18:2) moiety. In the ceramide subclass, Cer (d18:1/26:0) was connected to Cer (d18:1/26:1), as they both contained a sphingosine (d18:1) moiety. Clustering also depended on the nature of fatty acyl chain as exemplified by Cer (d18:0/16:0) which was connected to Cer (d18:1/16:0), as they both contained a palmitoyl moiety but displayed no connection with Cer (d18:1/18:0). To simplify the data visualization in Figure 6, the MN was organized on two orthogonal axes: the abscissa and the ordinate correspond to the number of insaturation of the sphinganine base moiety and to the fatty acyl chain length, respectively.

#### *2.6. Use of Retention Time Prediction for Lipid Annotation*

Applied to HCE cell lipidome, the MN approach was helpful to annotate more than 150 phospholipid and sphingolipid species (Table S3). Annotation was based on tandem mass spectrometry and confirmed by retention time with a maximum tolerance of 5% between theoretical and experimental tR values. In some cases, a co-elution and a co-selection of precursor ions in MS/MS was encountered leading to difficult mass spectra interpretation which did not readily permit to decide between two different lipid structures. In such a case, the tR prediction made it possible to annotate unequivocal lipid species.

For instance, the MS/MS spectrum of the ion at *m*/*z* 800.619 displayed the characteristic fragment ions at *m*/*z* 140.0123 and 196.0369 of a phosphoethanolamine headgroup suggesting a PE (40:1) structure. The spectrum also exhibited intense peaks corresponding to oleate (*m*/*z* 281.2480) and behenate (*m*/*z* 339.3254) but also two small peaks at *m*/*z* 253.2175 and 367.3488 indicative of palmitoleate and

lignocerate, respectively (Figure 6A,B). The selected precursor ion thus corresponded to a mixture of an abundant PE (22:0/18:1) and a less abundant PE (24:0/16:1). However, the presence of an ion at *m*/*z* 168.0444 could correspond to the headgroup of an isobaric [M−CH3] − precursor ion from PC (16:1/22:0). Thanks to the tR prediction, PC (16:1/22:0) was excluded as the difference between experimental and theoretical tR was 7% (Figure 6C). The origin of the phosphocholine ion at *m*/*z* 168.0444 was explained by the co-elution and co-selection by the Q1 quadrupole of the demethylated SM (d42:1) at *m*/*z* 799.668, a very scarce species in the precursor ion beam. This annotation was confirmed by a tiny difference of experimental and theoretical tR (0.3%).

**Figure 6.** Retention time prediction use to discriminate two isobaric PE and PC species. (**A**) MS/MS spectrum of PE (40:1) [M−H]<sup>−</sup> ion (*m*/*z* 800,619). (**B**) Low, middle, and high mass region of the spectrum A. (**C**) Theoretical tR plotted against experimental tR for standard PC and PE species (circles) and for PE (22:0/18:1) and PC (16:1/22:0) (squares).

#### *2.7. Use of Existing Lipid Library Database*

Although the proposed annotation procedure is performed through an individual inspection of each MS/MS spectrum, an automatization of the annotation may, at least in part, be considered using currently available lipid library especially LipidBlast or Lipidex [41–43].

For example, using LipidBlast, we performed annotation of the commercial standard lipids available in this study. Among the 65 lipids species studied, 50 were successfully annotated using LipidBlast and one misannotation was reported (Figure S4). None of the six SM species were annotated as SM subclass is not supported by the LipidBlast library using negative ion mode LC-ESI-MS/MS data. Moreover, among the 13 PC species included in the commercial standard lipid mixture, three PC failed to be annotated using LipidBlast.

In the frame of a lipidomic analysis, MN in combination with lipid databases may be regarded as a valuable and saving time approach giving a strong insight in the structural elucidation of lipid. Nevertheless, we believe that it remains important to use at least several known standard lipids to perform annotation of unknown lipid species with a high degree of confidence.

#### *2.8. E*ff*ect of HO on HCE Cells*

The lipid annotation through molecular networking and the retention time prediction approach proposed in this study was applied to assess lipid perturbations in human corneal epithelial cells exposed to hyperosmolarity (HO)—an *in vitro* model of dry eye disease [27,28]. Dry eye disease is a chronic inflammatory pathology of the ocular surface. It is characterized by alteration of tear film, HO, cell damage, and inflammation of the ocular surface, all contributing to a vicious circle [30,31,44,45]. Therapeutic strategies targeting inflammatory processes, such as cyclosporine or other anti-inflammatory agents, have been proposed to break this deleterious cycle [30,46]. To better understand the mechanism underlying DED and to find new marker of this pathology, especially to improve patient monitoring and to develop new targeted treatments, further investigations are still needed.

Hyperosmolarity is a key feature of DED [31]. Indeed, it induces significant stress targeting the corneal cell membrane [28,47,48]. Studying the modification of lipid homeostasis related to HO is relevant as this stressor induces the disruption of processes closely associated to cell membranes. Indeed, HO activates pro-inflammatory and pro-apoptotic processes, initiated at the cell membrane level. Hyperosmolarity stimulates downstream signaling pathways mediated by multiple membrane-bound proteins and enzymes [27,49–51]. The interest to study cell lipid disruption is also reinforced as HO favors ROS production which, in turn, targets lipids through peroxidation [52–54].

Lipidomic analysis associated to molecular networking and tR calculation led to the annotation of 150 lipid species and revealed that among them, 54 phospholipids and sphingolipids were significantly up- or downregulated in human corneal epithelial cells exposed to HO (Figure 7).

Regarding phospholipids, an increase of the cell concentration was observed for 4 PC, 10 PE, 3 PS, and 3 PI species under HO exposure. These phospholipids mainly contained oleate (18:1) located in the sn2 position. On the contrary, abundances of six ether-phospholipid species were strikingly decreased. Thanks to MN, they were identified as 2 PC-P, 2 PI-O, and 2 PE-P species containing polyunsaturated fatty acids at position sn2, especially arachidonic (20:4) and docosahexaenoic (22:6) acids (Figure 7). Ether phospholipids are known to be pools of FA (22:6) and FA (20:4) [55]. The FA (20:4) is the preferential substrate of cyclooxygenase (COX) and lipoxygenase (LOX), enzymes leading to pro-inflammatory eicosanoids. Therefore, our results suggest that HO, a key feature in DED, may favor the release of arachidonic acid from ether phospholipids to promote inflammatory process. Ether phospholipid, especially PE-P and PC-P species, are also known to be targets of oxidative stress through vinyl ether bonds. Beside inflammation, oxidative stress is also a key feature of DED pathophysiology [44,56]. Therefore, the significant decrease in PC-P and PE-P species observed in HCE cells under HO exposure may proceed through oxidation of these lipids which is known to generate toxic malondialdehyde (MDA) and 4-hydroxynonenal (4-HNE). Of note, MDA and 4-HNE release have previously been described in DED *in vitro* models as well as in conjunctival imprints of DED suffering patients [53,57].

Regarding sphingolipids, an increase of ceramide abundance was observed in HCE cells submitted to HO (Figure 7). Ceramides are bioactive lipid species known to promote inflammation through IL-1β release and to induce apoptosis through caspase 3 activation [7,58]. Dysregulation of ceramide metabolism is involved in many inflammatory diseases such as atherosclerosis, inflammatory bowel disease or multiple sclerosis [2,3]. In the frame of DED, apoptosis and inflammatory processes induced by ceramides may thus be considered as important mediators of the deleterious effects of HO in accordance with previous published data [47,59]. Thanks to MN, Cer (d18:0/16:0) and Cer (d18:1/16:0) were successfully identified. Under HO exposure, these two lipids species are increased in HCE cells. Because Cer (d18:0/16:0) and Cer (d18:1/16:0) are respectively substrate and a product of dehydroceramide desaturase—a key enzyme in *de novo* synthesis of ceramide—this result suggests that HO promotes *de novo* synthesis of ceramides making this ceramide biosynthetic pathway a putative therapeutic target in the frame of DED.

**Figure 7.** MN of lipids up- or downregulated in human corneal epithelial cells exposed to HO. MN displays (**A**) phospholipids and (**B**) sphingolipids species up- or downregulated in human corneal epithelial cells exposed to HO. The pie chart represents mean ion intensity of a lipid molecule from the control cell (yellow) and HO exposed cell (red) lipidomes. Of note, networking, based on the structural characteristics of lipid species, does not provide information regarding biosynthesis and modeling/remodeling pathways of the lipids found within HCE cells, especially under the effect of hyperosmolarity.

#### **3. Materials and Methods**

#### *3.1. Chemicals and Reagents*

Chloroform (Carlo Erba Reactifs SDS, Val-de-Reuil, France), acetonitrile, methanol, isopropanol, water of LC-MS grade (J.T. Baker, Phillipsburg, NJ, USA) and 3,5-di-*tert*-4-butylhydroxytoluene (Sigma–Aldrich, Saint-Quentin Fallavier, France) were used to perform cell lipid extraction and to prepare mobile phase for liquid chromatography. All commercial lipid standards were purchased from Avanti Polar Lipids, Inc. (Alabaster, AL, USA) and are listed in Table S1 of the Supplementary Materials.

#### *3.2. Sample Preparation*

The HCE cells were exposed to HO (500 mOsM) for 24 h and then were washed with Dulbecco's Phosphate-Buffered Saline (DPBS). Cells were harvested using trypsin-EDTA 0.05%, washed with DPBS, centrifuged at 2000 rpm for 10 min. Dry cell pellets were adjusted to 3 million cells and stored at −80 ◦C until analysis. After thawing, the cell pellets were resuspended in ultra-pure water (1 mL) and were sonicated for 5 min. Lipids were extracted using a chloroform/methanol/water (5:5:2, *v*/*v*/*v*) mixture containing 3,5-di-*tert*-4-butylhydroxytoluene 0.01% (*w*/*v*) as an antioxidant agent. Samples were subsequently centrifuged at 3000 rpm for 10 min, organic phases were collected, and solvents were evaporated under reduced pressure at 45 ◦C. Dry residues were dissolved in a 100 μL mixture containing acetonitrile/isopropanol/chloroform/water (35:35:20:10, *v*/*v*/*v*/*v*) before injection into the UHPLC-MS system.

#### *3.3. Data-Dependent LC-ESI-HRMS*/*MS Analysis*

Liquid chromatography-negative electrospray ionization mass spectrometry analysis of lipid extracts was performed on a UHPLC system (Waters®, Manchester, UK) combined with a Synapt®G2 High Definition MS™ (Manchester, UK) (Q-TOF) mass spectrometer (Waters®). Chromatographic separation was achieved on an Acquity® (Manchester, UK) CSH C18 column (100 <sup>×</sup> 2.1 mm; 1.7 μm). Lipids were eluted using a binary gradient system consisting in 10 mM ammonium acetate in acetonitrile/water mixture (40:60, *v*/*v*) as solvent A and 10 mM ammonium acetate in acetonitrile/isopropanol mixture (10:90, *v*/*v*) as solvent B. The eluent increased from 40% B to 100% B in 10 min, was held at 100% B for 2 min before returning to 40% B. The flow rate was kept at 0.4 mL.min<sup>−</sup>1. The column oven was set at 50 ◦C and the injection volume was 5 μL. The source parameters were as follows: capillary voltage 2400 V, cone voltage 45 V, source temperature 120 ◦C, desolvation temperature 550 ◦C, cone gas flow 20 L h<sup>−</sup>1, and desolvation gas flow 1000 L h−1. Leucine enkephalin (2 ng mL−1) was used as an external reference compound (Lock-Spray™, Manchester, UK) for mass correction. In a data-dependent acquisition mode (DDA), MS full scans were followed by MS/MS scans performed on the five most intense ions above an absolute threshold of 1000 counts. Selected parent ions were fragmented at collision energy ramp 20–40 eV and a selection window size of 1.0 Th. Scan durations for both MS and MS/MS were 0.2 s. In the full scan mode, the data were acquired between *m*/*z* 50 and 1200 using a resolution of 20,000 FWHM at *m*/*z* 500. Data acquisition was managed using Waters MassLynx™ software (version 4.1; Waters MS Technologies, Manchester, UK). A mixture of 65 standard lipids belonging to 9 of the main lipid classes (*i.e.*, phosphatidic acid (PA), phosphatidylethanolamine (PE), phosphatidylserine (PS), phosphatidylcholine (PC), phosphatidylglycerol (PG), ceramide (Cer), sphingomyelin (SM), HexosylCeramide (HexCer)) at a final individual concentration of 1 μM was also periodically injected throughout the analytical batch.

#### *3.4. Data-Preprocessing Parameters*

Raw data files were converted into universal open source mzXML file with MSConvert 3.0 and were then processed using MZmine 2.51 software. The MS and MS/MS spectra were extracted using MZmine 2.51 with a mass detection noise level set at 2E2 and 0E0, respectively. Chromatograms were then built with the ADAP algorithms [60] using a minimum group size of 5 scans, a group intensity threshold of 5000, and an *m*/*z* tolerance of 0.005 Da (about 10 ppm). The ADAP wavelets chromatogram deconvolution algorithm was used with the following settings: signal-to-noise ratio = 10, coefficient/area ratio = 50, peak duration range = 0.05−0.4 min, retention time wavelet range = 0.02–0.1, *m*/*z* range for MS/MS scan pairing of 0.01, and tR range for MS/MS scan pairing of 0.15 min. Chromatograms were de-isotoped using the isotopic peaks grouper algorithm with a *m*/*z* tolerance set at 0.005 (*m*/*z* < 500) and 10 ppm (*m*/*z* > 500), and a tR tolerance of 0.1 min. Peak alignment was performed using the join aligner method using the following parameters: *m*/*z* tolerance at 0.005 (*m*/*z* < 500) and 10 ppm (*m*/*z* > 500) and an absolute tR tolerance 0.15 min. Each MS/MS scans were associated with the corresponding MS scans using a tR tolerance of 0.1 min and a *m*/*z* tolerance of 0.005 (*m*/*z* < 500) and 10 ppm (*m*/*z* > 500). The peak list was finally gap-filled using the so-called module "same RT and *m*/*z* range gap filler" with *m*/*z* tolerance 0.005 (*m*/*z* < 500) and 5 ppm (*m*/*z* > 500).

#### *3.5. Molecular Network Analysis*

The MNs were created using the feature based molecular networking workflow of the Global Natural Products Social (GNPS) platform [61]. The following settings were used to build the network: minimum pairs Cos > 0.60, parent ion mass tolerance = 0.02 Da, fragment ion mass tolerance = 0.02, network topK < 100, minimum matched peaks = 6, and minimum cluster size = 2. The library spectra inquiries were performed using the same parameter values as those define for the network building. The MNs were finally visualized and annotated using Cytoscape 3.4.0 software (San Diego, California, USA) [62].

#### *3.6. Lipid Structure Assignment*

The structural annotation of unknown lipid species was based on the MNs generated on the GNPS platform using MS and MS/MS data as follow: (*i*) nodes associated to lipid standards were indexed using MzMine 2 thanks to MS, MS/MS data, and tR value, (*ii*) nodes associated to unknown lipids were subsequently indexed based on MS data, using online data base LIPIDMAPS and METLIN, and MS/MS data, through manual inspection on MzMine 2, (*iii*) annotation was finally supported by tR values by comparison of experimental tR values to the calculated one using tR prediction models. Based on these three criteria, lipids already annotated were used to create a database valuable for later annotation.

In order to demonstrate the relevance of the MNs in lipid species annotation, MS/MS spectra were individually inspected to select diagnostic product ions essential for annotation and the differences between theoretical and experimental *m*/*z* values were calculated using Excel software (see Supplementary Materials, compilation of experimental and theoretical *m*/*z* values of diagnostic product ions for PL and SL).

In accordance with the guidelines provided by the minimum reporting standards of the Metabolomics Standards Initiative [16], Table S3 includes the level of identification for all annotated lipids. Indeed, thanks to accurate *m*/*z* measurement, the MS/MS data inspection and retention time analysis, lipids annotated in HCE cells were assigned to group 1 or 2. Lipids for which we had the corresponding commercial standards were assigned to group 1. In addition, lipids for which we did not have the corresponding standards, the annotation was performed on the adequacy of *m*/*z* value, MS/MS data and retention time analysis and were thus assigned to group 2.

#### *3.7. Statistical Analysis*

A false discovery rate ((FDR)-adjusted *p* < 0.01) controlling procedure was performed to assess the statistical significance of the concentration differences of the identified lipids from cell extracts of HO-treated cells versus control cells. Each experiment was performed independently at least five times. The ANOVA, Dunnett test, and Student *t*-test were performed using GraphPad Prism 8 software (version 8; GraphPad Software, La Jolla, CA, USA) with a risk set at 0.05 (\* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001).

#### **4. Conclusions**

For the first time, the present study showed that the use of molecular networks makes it possible to facilitate and increase the reliability of lipid annotation in the course of lipidomic analysis. This approach, based on the use of the tandem mass spectrometry DDA mode in negative ionization conditions, allowed characterizing the fatty acyl chains of phospholipids and sphingolipids. In addition, if an ambiguity in the annotation of a lipid persists, the prediction of the retention time makes it possible to remove the latter. This new strategy makes it possible to cover the entire lipidome despite the limited number of standard lipids to which it is possible to have access commercially. The present study was limited to lipid subclasses which had structural characteristics which were clearly depicted by the MN under negative ionization conditions. Nevertheless, through this approach, we were able, in the context of a differential lipidomic analysis of an *in vitro* model of DED (*i.e.*, HCE cells exposed to HO) to annotate many lipids potentially involved in cell death and inflammation. Regarding others lipid subclasses such as glycerolipids, the positive ionization mode appears suitable to highlight their structural characteristics using molecular networks. It is currently in progress and will be the subject of a separate study.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2218-1989/10/6/225/s1, Figure S1: Extracted ion chromatograms of m/z values corresponding to precursor ions of standard phospholipids containing two palmitoyl moieties, Figure S2: MN of PS subclass, Figure S3: MN in PI subclass, Figure S4: Annotation of MS/MS spectra of commercial standard lipids using LipidBlast library, Table S1: Lipid composition of standard mix, Table S2: Repeatability and method precision (within-day, between-day, and intermediate precision), Table S3: Annotation of lipid species by MS/MS experiment. Supplementary Materials include the Excel table used for tR prediction and for diagnostic ion checking.

**Author Contributions:** Conceptualization, N.A., F.B.B., C.B., O.L., G.G.J., R.M., K.K., A.R., S.M.P.; methodology, R.M., N.A., G.G.J., F.B.B., K.K., A.R.; validation, N.A., G.G.J., F.B.B., K.K., A.R., S.M.P., O.L., C.B., R.M.; formal analysis, N.A., G.G.J., F.B.B., R.M., K.K., A.R.; writing—original draft preparation, R.M., N.A., F.B.B., K.K., A.R.; supervision, O.L., C.B., S.M.P., N.A.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by Sorbonne Université, the Institut National de la Santé et de la Recherche Médicale and Centre National de la Recherche Scientifique.

**Acknowledgments:** This work was completed with the support of the Programme Investissements d'Avenir IHU FOReSIGHT (ANR-18-IAHU-01). The authors thank the core facilities of the Institut de la vision, the Centre Hospitalier National d'Ophtalmologie des Quinze-Vingts, Région Ile-de-France and Ville de Paris.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Oxylipin Profiles in Plasma of Patients with Wilson's Disease**

#### **Nadezhda V. Azbukina 1, Alexander V. Lopachev 2, Dmitry V. Chistyakov 3,\*, Sergei V. Goriainov 4, Alina A. Astakhova 3, Vsevolod V. Poleshuk 5, Rogneda B. Kazanskaya 6, Tatiana N. Fedorova 2,\* and Marina G. Sergeeva 3,\***


Received: 17 April 2020; Accepted: 25 May 2020; Published: 29 May 2020

**Abstract:** Wilson's disease (WD) is a rare autosomal recessive metabolic disorder resulting from mutations in the copper-transporting, P-type ATPase gene ATP7B gene, but influences of epigenetics, environment, age, and sex-related factors on the WD phenotype complicate diagnosis and clinical manifestations. Oxylipins, derivatives of omega-3, and omega-6 polyunsaturated fatty acids (PUFAs) are signaling mediators that are deeply involved in innate immunity responses; the regulation of inflammatory responses, including acute and chronic inflammation; and other disturbances related to any system diseases. Therefore, oxylipin profile tests are attractive for the diagnosis of WD. With UPLC-MS/MS lipidomics analysis, we detected 43 oxylipins in the plasma profiles of 39 patients with various clinical manifestations of WD compared with 16 healthy controls (HCs). Analyzing the similarity matrix of oxylipin profiles allowed us to cluster patients into three groups. Analysis of the data by VolcanoPlot and partial least square discriminant analysis (PLS-DA) showed that eight oxylipins and lipids stand for the variance between WD and HCs: eicosapentaenoic acid EPA, oleoylethanolamide OEA, octadecadienoic acids 9-HODE, 9-KODE, 12-hydroxyheptadecatrenoic acid 12-HHT, prostaglandins PGD2, PGE2, and 14,15-dihydroxyeicosatrienoic acids 14,15-DHET. The compounds indicate the involvement of oxidative stress damage, inflammatory processes, and peroxisome proliferator-activated receptor (PPAR) signaling pathways in this disease. The data reveal novel possible therapeutic targets and intervention strategies for treating WD.

**Keywords:** COX; CYP450; LOX; oxylipins; PUFAs; lipidomics; UPLC-MS/MS; copper; Wilson's disease

#### **1. Introduction**

Wilson's disease (WD) is a rare autosomal recessive metabolic disorder resulting from mutations in the copper-transporting, P-type ATPase gene, ATP7B gene, which encodes a copper-transporting P-type ATPase [1]. The enzyme is responsible for the transport of copper into bile from hepatocytes, facilitating its incorporation into apoceruloplasmin to form ceruloplasmin, a major copper-transporting protein in the blood. The mutations lead to copper accumulation in the affected tissues [1–3], which

causes biochemical deviations, followed by the fluctuating occurrence of hepatic and extrapyramidal symptoms (EPSs), accompanied by the impairment of other organs (for more detail, see recent reviews [4–6]). The worldwide prevalence of WD is 1 in 30,000; it is even higher in populations with a high frequency of consanguinity [5]. The manifestations of WD are variable, and, in addition to liver diseases, may include neurological and/or psychiatric symptoms, as well as abnormalities in the blood or kidneys. It is therefore hypothesized that other genetic and/or environmental factors could influence the phenotypes of WD [5,7]. Besides different ATP7B mutations [8], other genetic variations may influence the variability in WD manifestation, such as apolipoprotein E (APOE), human prion protein (PRNP), 5,10-methylenetetrahydrofolate reductase (MTHFR), the interleukin-1 receptor antagonist (IL1RN), peroxisomal catalase, and other genes [7,9]. The influence of epigenetics, environment, age, and sex-related factors on the WD phenotype further complicates diagnosis [3,7,10]. The clear–phenotype correlations for WD are still unclear, although recent epigenetic whole genome screening data may shed light on in-depth pathogenic mechanisms [11]. Obtained from liver and blood samples from patients with WD, the data have shown specific sets of modified genes, enriched for functions in lipid metabolism and inflammatory responses [11]. Such changes can manifest themselves on the level of the organism as a whole in a variety of ways and be the cause of the observed differences in manifestations of the disease. This attracted attention to studying the disease on the level of the metabolome. Understanding the variations in the metabolome may help in identifying variations in the biochemical pathways leading to different manifestations of the disease. High-throughput techniques, such as metabolomic profiling, can deepen our understanding of the disease's pathogenesis and biology of WD manifestation [12,13], and therefore lead to new therapeutic approaches.

A promising type of metabolomic profiling is oxylipin measurement in the blood or other tissues. Oxylipins, derivatives of omega-3, and omega-6 polyunsaturated fatty acids (PUFAs) are signaling mediators that are deeply involved in innate immunity responses; the regulation of inflammatory responses, including acute and chronic inflammation; and other disturbances related to any system disease [14–17]. The conversion of PUFAs into oxylipins occurs via three major pathways, named according to their respective key pathway enzymes, such as the cyclooxygenase (COX), lipoxygenase (LOX), and cytochrome P450 monooxygenase (CYP450) branches of metabolism. Besides this, there are non-enzymatic conversions of PUFAs [15]. Due to the diversity of the individual oxylipin functions, it is difficult to predict the general direction of their action. For example, eicosapentaenoic (EPA) and docosahexaenoic (DHA) omega-3 PUFAs, as well as their derivative oxylipins, hydroxyeicosapentaenoic acids (HEPEs) and hydroxydocosahexaenoic acids (HDoHEs), are regarded as anti-inflammatory mediators [14,15]. Arachidonic acid (AA), an omega-6 PUFA, is mainly the source of prostaglandins (PGs), thromboxane (TX), leukotrienes (LTs), and hydroxyeicosatetraenoic acids (HETEs), attributed to groups of proinflammatory oxylipins. Meanwhile, cyclopentenone PGs, non-enzymatic metabolites of PGE2 and PGD2, possess anti-inflammatory features [18]. Oxidative derivatives of α-linolenic acid (ALA) can be transformed into hydroxyoctadecatrienoic (HOTrEs) acids or others [15]. Linoleic acid (LA)-derived oxylipins, such as hydroxyoctadecadienoic (HODEs) acids, agonists of PPARγ [19], or dihydroxyoctadecamonoenoic (DiHOMEs) acids, which are cytotoxic [20], exhibit both pro- and anti-inflammatory features [15]. Taken together, these data show that oxylipin synthesis should not be studied in groups of separate substances, but in terms of oxylipin profiles, which can characterize the different states of the studied organisms.

Indirect data indicate the possibility of oxylipin profile changes in WD. The roles of oxidative stress in the pathogenesis of WD, and dietary omega-3 PUFAs' usefulness in an animal model of WD suggest that oxylipins may be involved in the clinical manifestation of WD. Moreover, the number of some lipid-related nuclear receptors, such as retinoid X receptor (RXR), peroxisome proliferator-activated receptor α (PPARα), and hepatocyte nuclear factor 4 alpha (HNF4A), is generally decreased in WD animal models and humans [9,21–25]. It is important to note that although oxylipins exhibit multiple effects, they somehow change the state of the innate immunity system [13–15]. Oxylipins are important markers of activation of the system, including the regulation of inflammatory resolution processes [13,16]. Despite these facts, the role of the innate immunity system, and oxylipins as parts of the system, is still underestimated in many diseases.

The development of mass spectrometric methods for oxylipin detection made it possible to obtain oxylipin profiles from the plasma of patients with diseases, such as Alzheimer's disease [24], alcohol-related liver disease [25], or cancer [26]. However, no such studies have been conducted to characterize WD. Therefore, in the present study, UPLC-MS/MS lipidomics analyses were performed to characterize the plasma profiles of patients with various clinical manifestations of WD, compared to healthy subjects in order to identify the oxylipin characteristics of this disease.

#### **2. Results**

#### *2.1. Clinical Characteristics*

The study involved 39 WD patients and 16 healthy controls. The anthropometric, demographic, and blood biochemical parameters of the enrolled individuals are presented in Table 1. In total, 25 patients had the akinetic-rigid form, 10 had the trembling form, and 4 had other forms. Biochemical profiling of the blood of WS patients was conducted; data for ceruloplasmin and serum Cu concentrations are shown in Table 1.

**Table 1.** Demographic parameters of the patients, disease characteristics, and medications.


#### *2.2. Metabolomic Profiling*

Using UPLC-MS/MS, we detected a total of 43 metabolites in human plasma (Table S1). Metabolites were from different lipid classes: 3 PUFA (AA, DHA and EPA), 19 AA derivatives, one DGLA derivate, 7 DHA derivatives, 3 EPA derivatives, 7 LA derivatives, and 3 non-PUFA-derived compounds (OEA, AEA, Lyso-PAF).

#### *2.3. Volcano Plot Analysis*

To evaluate the separate metabolites that differ among WD and HC groups, we performed pairwise comparisons of age and gender-adjusted metabolite concentrations. The results were then illustrated using a volcano plot with Holm–Bonferroni correction (Figure 1). The four metabolites whose concentrations were changed significantly are indicated in red (12-HHT, EPA, PGE2, and PGD2). Barplots of the indicated compounds' relative concentrations are presented in Figure 1B.

**Figure 1.** (**A**) Volcano plot indicating significantly changed compounds. The X-axis indicates a log2 fold change of wilson disease WD to HC (healthy control) patients. Y-axis indicates −log10 *p*-values (adjusted). The cut-off for *p*-values is indicated based on Bonferroni correction. Compounds that changed insignificantly are indicated in gray, compounds whose means changed in WD (relative to HCs) more than twofold or less than twofold but insignificantly are indicated in green. Red dots stand for compounds, which changed more than twofold and had a *p*-value (adjusted < 0.05). (**B**) Relative concentrations of separate metabolites that changed significantly in WD patients in comparison with HCs. Pairwise comparison of adjusted means was conducted taking into account the age and sex of patients. \* *p* < 0.05 (adjusted for multiple testing).

#### *2.4. PLS-DA Model*

For data analysis, we used normalized concentrations of metabolites (see Section 2.6). The presence of outliers was identified by performing principal component analysis (PCA) to prevent their effects on the model. Hotelling's T2 test indicated three outliers in the healthy control group. A total of 52 samples, which were placed inside a 95% confidence interval ellipse red bounds (Figure 2A), were used for further analyses. For testing whether WD (Wilson disease) and HC (healthy control) patients could be distinguished based on oxylipin concentrations, the partial least square discriminant analysis (PLS-DA) was performed. The model was evaluated via cross-validation based on the overall error, balanced error rate (BER), and area under curve (AUC) values (Figure S1, Table S2). The optimal number of components was three. Projections on the first two components are presented in Figure 2B, and on the first three components in Figure 2C. Studied groups were separated with a small overlap. For each metabolite, the VIP score was estimated (as described in Section 2.6). The value of this parameter addresses the explained variation between classes in each projection. A total of seven metabolites, including 12-HHT, EPA, 14,15-DHET, 9-HODE, OEA, PGE2, and 9-KODE, with VIP score values > 1.5 are shown in Table 2.

**Figure 2.** (**A**) The principal component analysis (PCA) performed to verify outliers. The 95% Hotelling T2 confidence interval is indicated as an ellipse. (**B**) The partial least square discriminant analysis (PLS-DA) model discriminating healthy control (HC) and Wilson disease patients (WD). The explained variance of each component is indicated in brackets on the corresponding axis. (**C**) PLS2-DA model represented in 3-D showing separation among the HC and WD patients.

**Table 2.** Variable importance in projection (VIP) scores are shown for 7 metabolites. A cutoff value of 1.5 is established for VIP selection.


\* volcano plot indicating significantly changed compounds, *p* < 0.05 (adjusted for multiple testing).

#### *2.5. Similarity Matrix*

Since oxylipins represent different branches of metabolic pathways [14], we decided to estimate possible interconnections among compounds by calculating the pairwise association matrix, using data obtained using PLS-DA. A clustered image map (CIM), based on a hierarchical clustering of both the rows and the columns, was built using the Euclidean distance and complete linkage clustering algorithm (Figure 3). In the figure, each entry of the matrix is colored according to the association between metabolite concentrations and illness status (X and Y-variables in the model). The red color indicates positive correlation, whereas yellow/green indicates a weaker correlation. Dendrograms are shown on the left side (for metabolites) and on top (for patients). Color bar A indicates whether the patient belongs to WD (black) or HC (red). Based on the dendrogram and illness status, we subdivided the subjects into four groups (bar on the top of heatmap). WD patients can

be subdivided into three groups: Not distinguished from healthy donors (mix), with enrichment of HdOHE and HETE compounds, and with DiHETE, DiHOME enrichment (Figure 3).

**Figure 3.** A clustered image map was generated using the Euclidean distance and the complete linkage clustering algorithm. On the figure, each entry of the matrix is colored according to its value; rows represent metabolites, columns represent subjects. Dendrograms are shown on the left side (for patients) and on top (for metabolites). Color bars on the bottom of the picture indicate: (**A**) whether the subjects belongs to the Wilson disease (WD) group or healthy control (HC) group; (**B**) sex distribution: male (M) or female (F); (**C**) nephropathy status: no, yes, HC or data not available (N.A.); (**D**) psychosomatic status: no, yes, HC or N.A.; (**E**) form of the disease: trembling, akinetic-rigid, extrapyramidal, beforeneurology, or HC.

We further annotated the enrichment of modules using patients' clinical characteristics. The color bars on the bottom of Figure 3 indicate the clinical and anthropometric annotation of the patients. In total, five bars are presented on the CIM, indicating:


However, all WD patients clustered in a group on the left side turned out to be females (Figure 3). We tested associations with the Cu serum concentration, the severity of motor system dysfunction

according to the Shvab scale, age, and the debut age of subjects. There was no clear clustering according to the mentioned parameters (Figure S5A–D).

To test whether there were any differences in separate metabolites between selected groups, we conducted analysis of covariance (ANCOVA) to compare the adjusted means between groups taking into account the variability of the age and sex of patients. To identify which groups were different, pairwise comparisons of the adjusted means with the following Bonferroni multiple testing correction were applied. In the WD1 module patients, 10-HDoHE, 11-HETE, 12-HEPE, 12-HETE, 13-HDoHE, 15-HETE, 16-HDoHE, 5-HETE, and OEA were significantly different from HC (Figure S4). In the WD2 module patients, 10-HDoHE, 9-HODE, and AA were significantly different from HC. The mix module was not different from HCs (Figure S4).

#### *2.6. Pathway Enrichment Analysis*

Differences in separate metabolism branches are often a specific trait of biological processes [14–16]. This is why after independent analysis of the compounds, we took a step forward and investigated oxylipins as groups. Concentrations of compounds were summed according to their acid precursors (AA, DHA, EPA, ALA, DGLA, EA, EPA) (Figure 4A) or via the metabolic pathways they were derived from (cyclooxygenase (COX), cytochrome P450 monooxygenase (CYP), lipoxygenase (LOX), or non-enzymatic reactive oxygen species (ROS) (Figure 4B). It should be mentioned that in both cases, only the derivatives were summed up; free acids were grouped into the "others" unit. The classification used was in accordance with [15]. Then, a similarity matrix was calculated, and a complete linkage algorithm was performed for the acid precursor matrix (Figure 4A). To simplify the analysis between acid precursor and enzyme pathways, the second CIM was plotted using the order of the corresponding row as in Figure 4A, and clustering was performed only in columns (Figure 4B).

**Figure 4.** A clustered image map was performed using the Euclidean distance and complete linkage clustering algorithm. All polyunsaturated fatty acids (PUFA) derivatives were summed up according to their (**A**) initial substrate of biochemical pathways or (**B**) pathways' enzyme origin. In the figure, each entry of the matrix is colored according to its value, rows represent subjects, columns represent metabolites. Dendrograms are shown on the top (for metabolites). The color bar on the left side of the picture indicates whether a subject is WD (black) or a HC (red). Abbreviations: AA: arachidonic acid; DHA: docosahexaenoic acid, EPA: eicosapentaenoic acid; EA: AEA and OEA; LA: linoleic acid; DGLA: dihomo-γ-linolenic acid; ALA: α-linolenic acid; CYP: cytochrome P450 monooxygenase; LOX: lipoxygenase; ROS: reactive oxygen species; COX: cyclooxygenase.

The size of the cluster described in Section 2.5 increased. It was enriched with AA, DHA, EPA, LA, and DGLA metabolites, which indicates respective changes of their concentrations. It should be noted that the correlation between changes in the amount of metabolites of the CYP and LOX pathways was greater than that between their changes and changes in COX metabolites. However, this clustering was not explained by the clinical and demographic parameters mentioned in Section 2.4 and the reason for this standing apart is still unclear. Patients could be subdivided into two clusters independently of the grouping strategy (acids or enzymes), which is similar to the results presented in Figure 3. The bottom cluster of patients was associated with an overall upregulated content of oxylipins. On the other hand, the second cluster of patients did not show such an association.

After analysis of the grouped CIM, we determined a subgroup of patients as group 1. Patients from that group had significantly different concentrations of LOX metabolites (speaking about enzymatic pathways) and significant changes in the concentrations of AA, DGLA, and DHA derivatives and free fatty acids (Figure S5).

#### **3. Discussion**

Although WD is an autosomal recessive metabolic disorder, it possesses uncertain phenotype–genotype correlations and variability in its clinical manifestations. Understanding the biology of WD pathogenesis and improving diagnostic methods is an urgent problem posed before the science community [27]. Recently developed methods make possible quantitative measurement and analysis of a large number of markers, which, taken together, make up profiles. Oxylipin profiles are unique in that they reflect the activity level of a variety of biochemical processes in the organism and participate in the regulation of various signaling cascades. They possess a characteristic "fingerprint" when certain changes occur, and reflect dynamic characteristics of the organism, which is why interest in oxylipin profiles continues to grow.

Oxylipin profiles were investigated for the study of mechanisms, and as diagnostic markers for diseases, such as Alzheimer's disease [24], female breast cancer [26], alcohol-related liver disease [25], atherosclerotic diseases [28], and coronary artery disease [29]. Importantly, every disease was characterized by a special set of oxylipins: 9-HODE [26], 20-HETE [25], 8-HETE, LTB4, 9-HODE and 13-HODE [28], 9-HETE, and F(2)-isoprostanes [29]. Although t-test statistics of the oxylipin profiles of 39 (22 female and 17 male) WD patients and 16 (11 female and 5 male) donors allowed four substances (12-HHT, EPA, PGE2, and PGD2) to be revealed, which differ in WD vs HD, it is still not possible to suggest these substances as diagnostic markers, and further investigations are required. However, importantly, our data add new information concerning the biology of WD pathogenesis.

Indeed, analysis of data by VolcanoPlot showed that the patient vs. healthy groups differed significantly across three lipids. PLS-DA analysis revealed five more lipids that explained the difference in oxylipin profiles among WD and HC. Among them, two acids (EPA, OEA); two metabolites of LA, 9-HODE and 9-KODE, which can be attributed to LOX or non-enzymatic branches of metabolism; three metabolites of AA (12-HHT, PGD2, PGE2), attributed to the COX branch; and one AA metabolite from CYP branches of metabolism (14,15-DHET) were included. Although oxylipins possess multiple effects, and the same compound can be traced through various processes [14,15], the following processes can be characterized by the respective compounds: Oxidative stress (9-HODE, 9-KODE, ОЕA, ЕРA), inflammatory markers (9-HODE, 9-KODE, PGE2, 12-ННТ, PGD2), and peroxisome proliferator-activated receptor (PPAR) agonists (9-HODE, 9-KODE, OEA, EPA, 14,15-DHET). The data allow us to make assumptions about the possible signaling pathways involved in this pathology.

The ability of free Cu ions to participate in the formation of reactive oxygen species (ROS) and induce cellular toxicity is known [30]. In the presence of reducing agents (e.g., the superoxide anion radical), Cu2+ can be reduced to Cu+, which catalyzes the formation of hydroxyl radicals from hydrogen peroxide via the Haber–Weiss and Fenton reaction [31]. Therefore, the role of oxidative stress in the pathogenesis of WD is currently under investigation, and peroxisome impairment is suggested to be involved in WD pathophysiology [9,30,31]. Oxidized LA metabolites (HODE/KODE) are traditionally classified as oxidative stress markers [32]. Note that 9-HODE was also marketed as the most upregulated oxylipin species in the plasma of breast cancer patients, which indicates the possibility of oxidative stress involvement in this disease [26]. An interesting finding of our work is an increase of OEA and EPA, substances that are known preventers of oxidative stress [33,34]. Although oxidative stress can be viewed as a common disruption in various pathologies, which may suggest similarities in the oxidized forms of lipids, the variations in the oxylipin profiles obtained in various diseases do not support this point of view [24–26,28,29]. Our data sheds light onto some appropriate compensatory mechanisms that may be involved in WD.

An increase in 12-HHT, PGE2, and PGD2 points to the involvement of inflammatory processes in the pathogenesis of WD. This is in accordance with the data obtained in the animal model of WD (the Long-Evans Cinnamon rats), which is characterized by an increase in COX expression, the main enzyme in the synthesis of these substances [35]. Interestingly, dietary omega-3 PUFAs suppress acute hepatitis, prolong the survival of rats, and seem to lead to a decrease in COX expression [35]. Our data, showing an increase in 12-HHT, PGE2, and PGD2, are in accordance with this observation, because it is known that these omega-6 derivatives are inflammatory markers, which may be decreased by supplementation with dietary omega-3 PUFAs [36].

Importantly, the elevation of some oxylipins may lead to the activation of PPARs [37]. Three subtypes of PPAR (PPARα, PPARβ, PPARγ) are active regulators at the lipid metabolism and inflammation crossroad [38]. Recent studies have shown that PPARα and PPARγ are associated with steatosis and impairment of the antioxidant system in the liver of WD patients [39]. Inconsistent with this finding, a transcriptome analysis of the liver in the mouse model of Wilson's disease under copper-transporting, P-type ATPase gene Atp7b knockout identified the PPAR signaling pathway as a high-copper-responsive target pathway [40].

It is noteworthy that while PPARγ increased, PPARα mRNA expression is decreased with increased severity of WD [41]. This may reveal the "PPAR triad" mechanism that was conjectured for other cells, which respond to an excess of various types of PPAR ligands [42]. There are few data concerning the mechanisms of different PPAR-type changes in the presence of their ligands' excess at the organism level. Our data suggest that the increased number of PPAR ligands in the blood of WD patients may be associated with some kind of regulatory compensatory mechanism associated with the PPARs system. Our data also single out this signaling pathway for further consideration as being involved in the clinical manifestations of WD. Indeed, among the investigated substances, PPARα agonists were determined by OEA [41], EPA [37], and 14,15-DHET [43]. 9-HODE is an endogenous activator and ligand of PPARγ [19]. In this context, the question remains regarding the role of PGD2 and PGE2. Besides action via specific G-protein-coupled receptors, these prostaglandins are converted in the course of inflammatory reactions into prostaglandins 15d-PGJ2 and PGA2, respectively. Compounds with anti-inflammatory properties are formed, which activate PPARα and PPARγ [18]. At present, the question remains open whether an increase in PPAR agonists in the blood plasma of patients with WD is a protective mechanism that weakens the severity of the clinical course of the disease or, conversely, aggravates the symptoms. We were not able to find data on the use of synthetic agonists of PPAR in WD models. It is likely that fibrates and thiazolidinediones may be used as potential therapeutic agents for WD. Further research is required to elucidate the molecular mechanisms by which PPAR agonists may exert their effects in WD pathogenesis.

Beside PPARs, we cannot exclude oxylipins' involvement in the regulation of the activity of other nuclear receptors, because it is known that lipid-related nuclear receptors change in WD patients or animal models [22,23]. Decreased binding of the nuclear receptors FXR, RXR, HNF4α, and LRH-1 to promoter response elements and decreased mRNA expression of nuclear receptor target genes [22] and dysregulation of LXR/RXR heterodimers [23] may also be compensated by increased lipid agonist concentrations. Our data are consistent with these results, but further research is required to understand the mechanisms.

An important finding of our work is the subdivision of the patients into groups relative to the oxylipin profiles. Patients from the group had significantly different concentrations of LOX metabolites and had significant changes in the concentration of AA, DGLA, and DHA derivatives and free fatty acids. It is not yet clear which parameter underlies this division, since we did not find any correlation with gender or age. The subdivision into groups can reflect the various states of the innate immunity system. Along with various cytokines, oxylipins are part of the innate immunity system, being proinflammatory substances, as well as mediators of resolution [13]. Innate immune traits are more affected by the environment [44]. Environmental factors, including diet, exercise, stress, and toxins, profoundly impact the phenotypes of diseases, and WD is among them [7]. It is worthwhile to assume that the observed separation of patients by the oxylipin profile reflects a phenotypic response to environmental factors and therefore the current state of the innate immune system. This aspect requires further investigation of the metabolites that we found to be characteristic of WD.

#### **4. Materials and Methods**

#### *4.1. Reagents*

The oxylipins standards were as follows: 6-keto PGF1α-d4 (cat.no. 315210), TXB2-d4 (cat.no. 319030), PGF2α-d4 (cat.no. 316010), PGE2-d4 (cat.no. 314010), PGD2-d4 (cat.no. 312010), leukotriene (LT) C4-d5 (cat.no. 10006198), LTB4-d4 (cat.no. 320110), 5(S)-HETE-d8 (cat.no. 334230), 12(S)-HETE-d8 (cat.no. 334570), 15(S)-HETE-d8 (cat.no. 334720), oleoyl ethanolamide-d4 (cat.no. 9000552), EPA-d5 (cat.no, 10005056), DHA-d5 (cat.no. 10005057), and AA-d8 (cat. No. 390010) (Cayman Chemical, Ann Arbor, MI, USA). An Oasis®PRIME HLB solid-phase lipid extraction cartridge (60 mg, 3 cc, cat.no. 186008056) was obtained from Waters, Eschborn, Germany.

#### *4.2. Population and Study Design*

This was an observational study with 55 recruited people: 39 patients with WD and 16 healthy controls. In total, 39 individuals with Wilson's disease admitted to the regular inpatient treatment in the Research Center of Neurology (Moscow, Russia) were recruited for the study. Inclusion criteria for the WD patients included the following clinical and laboratory signs of the disease: Debut of the illness in childhood, adolescence, or adulthood (most often up to 35 years old); combined brain and internal organ damage (liver cirrhosis, hepatolienal syndrome, portal hypertension, tubular nephritis, etc.); damage to the central nervous system in the form of extrapyramidal syndrome; and systemic discuprinosis with impaired copper-ligand metabolism.

WD exclusion criteria included the following: Disease manifestation after 35 years; autosomal dominant type of inheritance; the presence of anamnestic, clinical, or paraclinical signs of another disease that can cause similar symptoms; hallucinations not related to medication; the presence of dementia or signs of impaired cortical function (aphasia, apraxia, etc.); the slowing down of vertical saccades or vertical gaze paralysis; a positive history of inflammatory diseases; chronic diseases and metabolic disorders; treatment with nonsteroidal anti-inflammatory drugs (NSAIDs) or corticosteroids during the last month; and pregnancy or breast-feeding during the study visit.

In total, 16 healthy individuals not affected by neurodegenerative disorders as verified by clinical examination were included in the study. They were recruited among people undergoing periodic health examinations at the same center. The exclusion criteria for healthy controls were the same as for patients with WD.

The Ethics Committee of the Research Center of Neurology approved this study (protocol №4-4/19 15.05.19), and informed written consent was obtained from each patient and control according to the guidelines approved under this protocol (Article 20, Federal Law "Protection of Health Right of Citizens of Russian Federation" N323- FZ, 11.21.2011).

#### *4.3. Clinical Evaluation*

The criteria for inclusion in the group of patients with WD consisted of clinical and laboratory signs of the disease: Debut in childhood, adolescence, and adulthood (most often up to 35 years old); combined damage to the brain and internal organs (cirrhosis, hepatolienal syndrome, portal hypertension, tubular nephritis, etc.); central nervous system (CNS) damage in the form of extrapyramidal syndrome, including tremor, stiffness, dysarthria, dysphagia, cognitive impairment, and dysphoria; and extra-neural symptoms, including hepatosplenomegaly, hemorrhages, failure of levodopa treatment (prescribed in connection with Parkinson's syndrome). Criteria for laboratory diagnosis of WD: Kaiser–Fleischer corneal ring (when using a slit lamp); a decrease in the concentration of ceruloplasmin copper ligand protein in the blood serum; hypersecretion of copper with urine; increase in the concentration of free copper in blood serum; a decrease in the concentration of total copper in serum; decreased serum zinc concentration; increased copper concentration in liver biopsy specimens; DNA diagnostics (detection of mutations in the ATP7B gene); and high therapeutic effect when using copper eliminating chelates (D penicillamine, trientin) and zinc preparations. An additional diagnostic criterion for WD was the result of neuroimaging (computed tomography and magnetic resonance imaging (CT, MRI)), showing an atrophic process in the cerebral hemispheres, cerebellum, subcortical structures with a corresponding expansion of subarachnoid spaces, and the ventricular system, as well as foci in the area of lenticular nuclei, globus pallidus, and visual hillock [26,27].

#### *4.4. Blood Sample Collection*

Taking into account the reported serum oxylipin variety during daytime [45], all blood sample collection was conducted in the morning in the fasted state. The plasma was obtained immediately after blood sampling, aliquoted, and stored at −80 ◦C for further analysis.

#### *4.5. UPLC-MS*/*MS Conditions and Sample Preparation*

Samples were prepared for MS analysis by the solid-phase extraction (SPE) method using an Oasis®PRIME HLB cartridge (60 mg, 3 cc). For eicosanoid extraction, plasma (900 μL) was deproteinized with 1 mL of methanol, vortexed, and centrifuged at 12,000 rpm for 5 min at ambient temperature. The supernatant was diluted 1:6 with mQ water containing 0.1% formic acid for the next steps of SPE. Then, the sample was loaded, and the cartridge was washed with 2 mL of 15% methanol containing 0.1% formic acid, after which the lipids were sequentially eluted with 500 μL of anhydrous methanol and 500 μL of acetonitrile. The resulting samples were concentrated by evaporation of the solvent under a gentle stream of nitrogen and stored at −80 ◦C. For the identification of lipid mediators, the respective lipid extracts were analyzed using an 8040 series UPLC-MS/MS mass spectrometer (Shimadzu, Japan) in multiple-reaction monitoring mode at a unit mass resolution for both the precursor and product ions [46]. The selected molecular ions were fragmentized in the gas phase by collision-induced dissociation and analyzed by tandem (MS/MS) mass spectrometry. The studied metabolites were identified and quantified according to the comparison of their multiple reaction monitoring parameters, retention times, and peak areas with the parameters obtained for deuterated internal standard compounds of the same classes (6-keto PGF1α-d4, TXB2-d4, PGF2α-d4, PGE2-d4, PGD2-d4, leukotriene (LT) C4-d5, LTB4-d4, 5(S)-HETE-d8, 12(S)-HETE-d8, 15(S)-HETE-d8, oleoyl ethanolamide-d4, EPA-d5, DHA-d5, AA-d8) (Table S1) using a commercial software method package Lipid Mediator Version 2 (Shimadzu, Tokyo, Japan) according to the manufacturer's instructions.

Prior to analysis, the plasma oxylipin detection method was validated according to food and drug administration (FDA) recommendations [47]. The stock ethanol solution containing 2 ng of each of the 15 deuterated oxylipin standards was prepared. Calibration samples (1.4, 1, 0.4, 0.2, and 0 ng/probe) were obtained by further dilution of the stock solution containing the same total volume of plasma. For each standard intraday, the interday reproducibility and relative standard deviation (RSD, %) were determined. The accuracy was measured with spiked standards for 3 concentration ranges:

0.2–0.8 ng/probe, 0.9–1.3 ng/probe, and 1.4–2 ng/probe. The limit of detection (LOD) was determined as the signal to noise ratio = 3, and the limit of quantification (LOQ) as the signal to noise ratio = 3. Signal to noise ratio was calculated as the standard error of the regression/slope. Results are presented in Table S2.

#### *4.6. Experimental Data Analysis and Statistics*

Comparison of the relative concentrations was performed using the two-sample two-sided *t*-test, followed by Bonferroni–Holm correction for multiple comparisons. *p* < 0.05 was considered as statistically significant.

Metabolomics data was analyzed using the mixOmics R package version 6.1.1 [48]. After data normalization on internal standards, peak area mean centering and unit variance scaling was applied. Class separation was analyzed by partial least square discriminant analysis (PLS-DA). The quality of the built model was estimated using leave-one-out cross validation. Each round of cross-validation included training on the bigger data subset and validation on the randomly selected sample. The model's predictive performance was estimated based on the validation results combined over rounds. This procedure was repeated for different numbers of components in the PLS-DA model. The overall error, balanced error rate, and AUC were obtained for each number of components (Figures S1 and S2, Table S5).

After building the PLS-DA model with 3 components, VIP scores for each investigated metabolite were calculated. A VIP score is a weighted sum of squares of the PLS loadings regarding the explained variation in each projection. A cutoff for VIP-scores was accepted as 1.5 according to the metabolomics standard initiative (level MSI = 1).

Analysis of covariance (ANCOVA) was used to compare the means of single metabolites between studied groups, taking into account sex and age. ANCOVA was performed using rstatix package for R. Pairwise comparisons of relative metabolites' concentrations was performed using function emmeans\_test (rstatix package), also taking into consideration age and sex as covariates. Analysis was followed by Bonferroni–Holm correction for multiple comparisons. *p* < 0.05 was considered as statistically significant.

#### **5. Conclusions**

In conclusion, our findings reveal alterations of the plasma oxylipin profiles in Wilson's disease patients, and heterogeneity in patients relative to the oxylipin profiles. Eight lipids were found to vary between HC and WD: EPA, OEA, 9-HODE, 9-KODE, 12-HHT, PGD2, PGE2, and 14,15-DHET; among them, PGE2, PGD2, 12-HHT, and EPA changed significantly (based on the pairwise comparison of means adjusted for sex and age).

The biological significance of these compounds indicates the involvement of oxidative stress damage, inflammatory processes, and PPAR signaling pathways in this disease. The data reveal novel possible therapeutic targets and intervention strategies for treating WD.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2218-1989/10/6/222/s1, Figure S1: The Balanced Error Rate (BER) and overall error rate estimated via cross-validation for different component numbers, Figure S2: ROC (receiver operating characteristic) curve for chosen number of components, Figure S3: A clustered image map was generated using euclidean distance and complete linkage clustering algorithm, Figure S4: Relative concentrations of separate metabolites which changed significantly in WD1 and/or WD2 in comparison with HC, Figure S5: Relative concentrations of summed metabolites which changed significantly between group1 and group2, Table S1: UPLC-MS/MS parameters of the identified lipids, Table S2: Intraday reproducibility, Table S3: Interday reproducibility, Table S4: Accuracy, LOD, LOQ, Table S5: AUC (Area Under Curve) values and p-value estimated for different number of components, Table S6: Mean +/− standard deviation of relative concentrations for healthy controls (HC) and Wilson disease patients (WD), Table S7: Source acid and metabolic enzyme for the analyzed oxylipins.

**Author Contributions:** Experimental procedures, N.V.A., A.A.A., S.V.G., D.V.C., V.V.P. and A.V.L.; writing—original draft preparation, N.V.A., D.V.C. and M.G.S.; writing—review and editing, A.V.L., M.G.S., R.B.K and T.N.F.; clinical analysis, A.V.L., V.V.P., R.B.K. and T.N.F.; performed the MS analysis, D.V.C. and S.V.G.; data analysis, N.V.A.; supervision, T.N.F. and M.G.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The reported study was funded by RFBR according to the research project № 19-29-01243.

**Acknowledgments:** The publication has been prepared with the support of the "RUDN University Program 5-100" (MS/MS analysis).

**Conflicts of Interest:** The authors declare no competing interests.

#### **Abbreviations**


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Lipidomic Phenotyping Reveals Extensive Lipid Remodeling during Adipogenesis in Human Adipocytes**

**Florian Miehle 1,2, Gabriele Möller 1, Alexander Cecil 1, Jutta Lintelmann 1, Martin Wabitsch 3, Janina Tokarz 1, Jerzy Adamski 1,2,4,5 and Mark Haid 1,\***


Received: 28 April 2020; Accepted: 23 May 2020; Published: 26 May 2020

**Abstract:** Differentiation of preadipocytes into mature adipocytes is a highly complex cellular process. At lipidome level, the adipogenesis remains poorly characterized. To investigate the lipidomic changes during human adipogenesis, we used the LipidyzerTM assay, which quantified 743 lipid species from 11 classes. The undifferentiated human SGBS cell strain showed a heterogeneous lipid class composition with the most abundant classes, phosphatidylethanolamines (PE), phosphatidylcholines (PC), and sphingomyelins (SM). The differentiation process was accompanied by increased ceramide concentrations. After completion of differentiation around day 4, massive lipid remodeling occurred during maturation, characterized by substantial synthesis of diacylglycerols (DAG), lysophosphatidylethanolamines (LPE), PC, PE, SM, and triacylglycerols (TAG). Lipid species composition became more homogeneous during differentiation to highly concentrated saturated and monounsaturated long-chain fatty acids (LCFA), with the four most abundant being C16:0, C16:1, C18:0, and C18:1. Simultaneously, the amount of polyunsaturated and very long-chain fatty acids (VLCFA) markedly decreased. High negative correlation coefficients between PE and PC species containing VLCFA and TAG species as well as between ceramides and SM imply that PE, PC, and ceramides might have served as additional sources for TAG and SM synthesis, respectively. These results highlight the enormous remodeling at the lipid level over several lipid classes during adipogenesis.

**Keywords:** adipocytes; adipogenesis; differential mobility spectrometry (DMS); lipidomics; lipidyzer; mass spectrometry; metabolomics; phenotyping; Simpson-Golabi-Behmel syndrome (SGBS)

#### **1. Introduction**

Overweight and obesity have increased dramatically in recent decades and now affect hundreds of millions of people worldwide, reaching pandemic levels [1]. Obesity has an adverse impact on a range of physiological processes, thereby increasing the risk of developing diseases like type 2 diabetes [2], cardiovascular diseases [2,3], and some types of cancer [4]. Overweight and obesity are mostly characterized by an excess of white adipose tissue (WAT). Adipocytes, the main constituent of this tissue, control the energy balance by storing triacylglycerols in periods of energy excess and breaking down these lipids during energy deprivation. However, the physiological role of adipocytes is much more complex than simply acting in energy storage. These cells secrete numerous diverse lipids and proteins controlling and regulating various bodily functions like appetite, immunological and inflammatory responses, and blood pressure and thereby act as an endocrine organ [5,6].

The development of adipocytes from precursor cells is known as adipogenesis. Within this differentiation process, fibroblast-like preadipocytes differentiate into lipid-laden and insulin-responsive adipocytes. This highly complex process involves the concerted interaction of a cascade of transcription factors like peroxisome proliferator-activated receptor gamma (PPARγ) and CCAAT/enhancer-binding proteins (C/EBPs) as well as different metabolic pathways including the TCA cycle, fatty acid synthesis, glycolysis, and polyamine biosynthesis [7–9]. The differentiation process has been characterized predominantly in murine cells, using different omics approaches like transcriptomics [10,11], metabolomics [9,12], proteomics [13–15], as well as the combination of transcriptomics and metabolomics [16]. However, lipids and their precise composition changes were so far sparsely characterized due to a lack of appropriate "high-resolution" mass spectrometry (MS) methods and internal standards. Several years ago, Roberts and coworkers characterized the levels of free fatty acids as well as the total levels of fatty acids of triacylglycerols (TAG) and glycerophospholipids in differentiating murine 3T3-L1 cells [9]. Two other studies analyzed the levels of some phosphatidylcholines (PC) using the same cell model, but without resolving the lipid isobars in experimental setups that were focused on the analysis of polar analytes [16,17]. Liaw and coworkers compared differentiated 3T3-L1 cells with primary mouse ear-derived mesenchymal stem cells and brown BAT-C1 adipocytes with a global lipid profiling method [18]. Additionally, they investigated the differentiation of murine adipocytes. However, they did not longitudinally analyze the whole adipogenesis process but compared only undifferentiated with fully differentiated 3T3-L1 cells. Even less knowledge on lipids has been compiled in the past for human adipocytes. Collins and coworkers analyzed the levels of fatty acid compositions of TAG and phospholipids without further class separation in primary adipocytes of subcutaneous origin [19]. To the best of our knowledge, the adipogenic process in a human cell model was so far not characterized with high-resolution lipidomics approaches. Therefore, we studied the development of human preadipocytes into mature adipocytes on a lipidomics scale with the recently developed LipidyzerTM (SCIEX, Darmstadt, Germany) method. We were able to simultaneously quantify 743 lipid species from 11 lipid classes. In combination with multivariate statistics, we could uncover correlations that suggest that extensive lipid remodeling occurs between several lipid classes during adipogenesis. These findings might contribute to the elucidation of new therapy strategies in obesity and other lipid metabolism affected disorders.

#### **2. Results**

#### *2.1. Analytical Method Validation*

In order to study the process of adipogenesis of human SGBS cells we decided to use the novel targeted Lipidyzer™ technology. Although originally developed and validated for the fast and automated analysis of human plasma samples [20], recently published studies also showed good performance of the Lipidyzer™ assay with other matrices than plasma [21–25]. With the aim to obtain a meaningful "fit for purpose method", we investigated the analytical performance regarding linearity and repeatability for the SGBS cells.

For linearity evaluation, we used nine different volumes of cell homogenates (10, 20, 40, 50, 60, 80, 100, 200, and 300 μL) of undifferentiated and differentiated (day 15) SGBS cells and determined mean concentrations for the single lipid classes. For most classes, the coefficients of determination (R2) of the linear regressions were higher than 0.9 for both sample types, the differentiated and the undifferentiated cells, meaning we had good linearity of the method (Figure S1, Table S1). However, two classes, the dihydroceramides (DCER) and free fatty acids (FFA), revealed insufficient linearity for both sample types. The measurements of lactosylceramides (LCER) showed high linearity in the undifferentiated (*R*<sup>2</sup> = 0.9734) but low linearity in the differentiated cells (*R*<sup>2</sup> = 0.0479).

For repeatability evaluation, we calculated the relative standard deviations (CV) of the QC pooled samples. With the exception of DCER, all other lipid classes revealed CV values < 15% (Table S2). We also analyzed the background signals of the lipid classes by defining a lower cut-off value for the Lipidyzer™ method at a value of 1.5× the value measured in the blank. In total, 12 lipid classes were present in QC samples in quantities higher than 1.5× the blank values (Table S2). Only FFA showed up in lower apparent amounts because their measured concentrations were already high in blank samples. In conclusion, the Lipidyzer™ assay showed good linearity and repeatability for most of investigated lipids in SGBS cells.

Finally, we decided to exclude the FFA and DCER data from the data set due to their insufficient reliability in the validation testing. However, we left the LCER in the data set as we think that the drop of concentration between undifferentiated and differentiated cells is a noticeable result of our study.

#### *2.2. Cellular Lipid Composition Undergoes Remodeling During Adipogenesis to Mainly TAG*

Alterations in the lipid content of some lipid classes (e.g., TAG, phospholipids) in murine cells undergoing the process of adipogenesis have long been known [9,16,18,19]. However, concentration changes of many lipid species from multiple lipid classes at the different stages of human adipogenesis have not been investigated so far. We used the LipidyzerTM technology to follow the differentiation process of the human SGBS cell strain by quantifying the lipids of samples at days 0, 4, 8, 12, 16, and 20 of adipogenesis.

To track the successful cell differentiation of preadipocytes into lipid-laden adipocytes, we monitored the cellular process by microscopy (Figure S2) and analyzed relative mRNA levels of the main adipogenic transcription factors *PPARG* and *CEBPA* (Figure S3). Microscopic analysis showed enormous lipid storage in droplets starting between day 4 and 8. The analyses of the relative mRNA expression levels showed strong upregulation of *PPARG* (40.7 ± 9.5 fold change at day 12 compared to day 0) and *CEBPA* (52.0 ± 8.8 fold change at day 12 compared to day 0). These data demonstrate successful differentiation of SGBS cells into adipocytes.

We were able to simultaneously quantify 743 lipid species of 11 different lipid classes with the accurate identification of lipid isobars. To investigate putative differences between the lipid concentration levels at different time points, we conducted partial least squares-discriminant analysis (PLS-DA, Figure 1) and principal component analysis (PCA, Figure S4). PLS-DA shows a clear separation of the different time points of adipogenesis using the first two principal components with 68.5% and 12.7% of explained variance (Figure 1). While component 1 was sufficient to separate the early phase of SGBS differentiation (days 0, 4, and 8), the second component was necessary for separation of the later stages of differentiation (days 12, 16, and 20). Component 1 was mostly influenced by TAG species, whereas their influence on component 2 was lower (Table S3).

Next, we were interested in the lipid class compositions at the six time points of differentiation (Figure 2A). In preadipocytes, that is, at the start of differentiation (day 0), a very heterogeneous lipid composition was observed with the most dominant classes being the phosphatidylethanolamines (PE; 32.1% ± 0.9%), phosphatidylcholines (PC; 26.3% ± 0.4%), sphingomyelins (SM; 19.7% ± 0.6%), and TAG (10.3% ± 0.5%). While the relative proportions of PE, PC, and SM declined tremendously during the ongoing cell differentiation, the fraction of TAG increased from initially 10.3 ± 0.5% to finally 96.9 ± 0.4% at day 20. The relative fractions of all other lipid classes decreased strongly during the

differentiation process to fraction sizes of finally 1.2% and lower. To conclude, the different stages of cell differentiation could be clearly distinguished based on the relative lipid compositions (Figure 2A) as well as by PLS-DA (Figure 1). Interestingly, the relative lipid class compositions did not reveal strong changes after day 8 of differentiation.

**Figure 1.** Partial least squares-discriminant analysis (PLS-DA) score plot showing very clear clustering of lipid species regarding the different time points of adipogenesis. While component 1 (68.5% variance) was sufficient to separate the early phase of differentiation (days 0, 4, and 8), the second component was necessary for separation of the later stages of differentiation (days 12, 16, and 20, 12.7% variance). The color code for the data points of the different days of adipogenesis is shown in the box inside the figure. Illustrated are also the 95% confidence intervals for each group. Cross-validation and permutation results confirmed the model to be predictive and not overfitted (R2: 0.97, Q2: 0.97; *<sup>p</sup>* <sup>&</sup>lt; <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> (0/2000 permutation numbers), test statistics selected by separation distance (B/W)). Each group consisted of six samples.

**Figure 2.** Changes in lipid class compositions and concentrations in specific lipid classes in the course of adipogenesis showed enormous lipid remodeling during the ongoing differentiation process. (**A**): Relative lipid class compositions at different days of adipogenesis in molarity %. Prior to differentiation at day 0, a very heterogeneous class distribution could be observed. This composition became more uniform during adipogenesis due to the predominance of triacylglycerols (TAG). The relative fraction for TAG increased from 10.3 ± 0.5 molarity % at day 0 to 96.9 ± 0.4 molarity % at day 20. (**B**): The concentration profiles of different analyzed lipid classes showed individual time courses during ongoing adipogenesis. Cholesteryl esters (CE), HCER, and LCER concentrations decreased strongly during adipogenesis, whereas LPE, PC, PE, and TAG concentrations strongly increased. Ceramide (CER) concentrations increased from day 0, peaking at day 4, and then decreased below the starting concentration levels. DAG species had their concentration maximum from days 8 to 12. LPC fluctuated during the whole differentiation process.

The time courses of the lipid species concentrations revealed an interesting alternative perspective on adipogenesis (Figure 2 B). Overall, we observed significant changes over time in the concentrations of 725 out of 743 lipids, that is, for 97.7% of all quantified lipid species. When looking at the lipid classes, only three of them, namely, CE (*<sup>p</sup>* <sup>=</sup> 4.52 <sup>×</sup> <sup>10</sup><sup>−</sup>4, Table S4), HCER (*<sup>p</sup>* <sup>=</sup> 5.23 <sup>×</sup> <sup>10</sup><sup>−</sup>3), and LCER (*p* = 2.80 <sup>×</sup> 10<sup>−</sup>5), showed continuous and significant decreases in concentration levels from day 0 to day 20 (Figure 2B). CER concentrations increased significantly from day 0 to day 4 of differentiation (*<sup>p</sup>* <sup>=</sup> 1.51 <sup>×</sup> <sup>10</sup><sup>−</sup>5) and subsequently decreased significantly below the limit of detection (LOD). In contrast, the TAG concentration levels increased continuously (*<sup>p</sup>* <sup>=</sup> 4.01 <sup>×</sup> <sup>10</sup><sup>−</sup>6). The TAG concentration levels were the highest among all classes and levels started from 19.2 ± 2.0 μM on day 0 to 8874.3 ± 1072.1 μM on day 20. LPE, PC, and PE concentrations also increased strongly during adipogenesis, whereas SM concentrations showed only a mild increase during differentiation (*p* = 5.04 <sup>×</sup> 10−5). In contrast, DAG concentrations increased from day 4 to day 8 and remained at a high level (*p* = 3.40 <sup>×</sup> 10−5). LPC concentrations fluctuated around the starting value during differentiation (*<sup>p</sup>* <sup>=</sup> 1.50 <sup>×</sup> <sup>10</sup><sup>−</sup>4).

#### *2.3. The Most Abundant Fatty Acids in Di*ff*erentiated Human SGBS Cells Are C16:0, C16:1, C18:0, and C18:1*

As we were interested in the concentration changes of the single fatty acid (FA) species during SGBS adipogenesis, we subsequently focused on a detailed analysis of the fatty acids bound to the lipid backbone. Figure 3 illustrates the time courses of the concentrations for the single FA side chains summarized over all lipid classes. With the exception of C20:4, all LCFA as well as the medium-chain FA (MCFA) lauric acid (C12:0) showed strong increases in concentration during adipogenesis. Among this group of FA, C18:1, C16:0 (palmitic acid), C16:1, and C18:0 (stearic acid)—in descending order—were the most abundant. In contrast to the MCFA and nearly all LCFA, most of the VLCFA decreased during adipogenesis. FA C22:0, C22:1, and C22:2 showed fluctuating concentration profiles.

A detailed analysis of the time courses of the concentrations for the single FA side chains separately for each lipid class allowed an even more detailed glimpse into adipogenesis (Figure 4). Illustrations on an enlarged scale can be found in Figure S5A–F. We were able to identify individual FA concentration changes over time of differentiation, which were strongly dependent on the lipid class.

The FA species of the three classes of ceramides, namely CER, HCER, and LCER, showed similar behavior in that the concentrations of the LCFA and VLCFA side chains decreased during adipogenesis. In contrast, the LCFA in DAG, LPE, PC, PE, and TAG substantially increased. In particular, the concentrations of C18:1 and C16:0 increased strongly during adipogenesis; for example, TAG-C16:0 increased from 2.5 μM at day 0 to 1515.1 μM at day 20. In contrast, VLCFA were present at only very low levels in the classes of DAG, LPE, PC, PE, and TAG.

The FA composition of the SM differed strongly from the compositions of the other classes in that the SM comprised the highest absolute amounts of VLCFA. Especially C22:0 (behenic acid) and C24:0 (lignoceric acid) were present in SM at considerable concentrations. Their concentrations were increased up to a factor of 320 at day 20 compared with the other classes. The fatty acid concentrations of the LPC did not change markedly during adipogenesis. In contrast, the LCFA of the CE decreased to half the maximal concentrations and the VLCFA even more during the differentiation process.

**Figure 3.** Time courses for the concentrations of the single FA side chains bound in the lipids summarized over all classes. Analysis of summarized lipids' side chain concentrations over all classes revealed four dominant fatty acids, namely FA 18:1, FA 16:0, FA 16:1, and FA 18:0. Their concentrations increased strongly during adipogenesis together with the other LCFA (with the exception of FA 20:4) and the MCFA 12:0. In contrast, the very VLCFA mostly decreased during adipogenesis.

**Figure 4.** Fatty acid concentrations and compositions changed markedly during adipogenesis in all 11 lipid classes. At the start of differentiation, the lipids had a very heterogeneous side chain distribution with high concentration levels of LCFA and VLCFA. The concentrations of the VLCFA decreased during adipogenesis in all classes, with the exception of the class SM. The concentration courses of the LCFA were more complex because their levels increased markedly during adipogenesis within the classes DAG, LPE, PC, PE, and TAG, but decreased in the class CE. The FA concentration course of the SM differed strongly from the other classes because the VLCFA remained at high levels during adipogenesis.

#### *2.4. Correlations between Concentration Profiles of Lipid Species from Di*ff*erent Lipid Classes Reveal Extensive Lipid Remodeling during Adipogenesis*

The opposing trends in concentrations of some lipid classes raised the questionwhether they might be the result of an underlying regulatory network. To reveal possible associations between the lipid species, we computed pairwise Spearman's rank correlations of lipid species concentration trajectories (Figure 5). The lipids could be assigned to six clusters (Figure 5A) and for each of the clusters the average concentration changes over time are displayed in Figure 5B. Species of different lipid classes were found to be distributed between the different clusters (Figure 6A and Figure S6).

**Figure 5.** Spearman's rank correlation analysis of lipid concentration trajectories during adipogenesis showed strong clustering and correlation of the lipid species. Panel (**A**) illustrates the matrix of the analysis where each square indicates the Spearman's rank correlation coefficient. Positive correlations between the variables are shown in red, while negative correlations are shown in blue. Correlation matrix enables the assignment of six lipid clusters. Species were clustered using Ward's clustering algorithm. Panel (**B**) shows the changes of the average lipid concentration in the different clusters with time.

Cluster 1 (*n* = 97 lipid species) consisted of lipid species with continuously decreasing concentration profiles (Figure 5B). In cluster 1, lipid species with the highest average FA side chain lengths (expressed as average total number of C-atoms) as well as number of double bonds (DB) were found, independent of the lipid class (Figure 6B–D). The cluster was dominated by PE (37.6%) and PC (20.8%) species. However, with the exception of SM, lipids from all other lipid classes were also clustered here (Figure 6A). Interestingly, 80% of all LCER species of the dataset could be found in this cluster (Figure S6A). In addition, more than 74.3% of the lipid species within this cluster had at least one FA side chain with at least 20 C-atoms and a high degree of desaturation (Table S5). Therefore, this cluster can be considered as the "PUFA cluster".

Cluster 2 (*n* = 61) comprised species that exhibited decreasing concentrations starting at day 8. This cluster contained lipids of all classes except those of LPE. The species had approximately the third highest total number of C-atoms and DB (Figure 6B–D). Remarkably, cluster 2 contained a relatively high number of sphingolipids: 56.6% of all CER, 40.0% of all HCER, and 41.7% of all SM species of the dataset (Figure S6B) could be found in this cluster. Therefore, it can be characterized as the "sphingolipid and PUFA cluster".

Cluster 3 species (*n* = 19) showed a fluctuating concentration course with a decrease in concentrations until day 8, followed by an increase to above the starting values. Four lipid classes were part of this cluster: TAG (42.1%), DAG (10.5%), PC (15.8%), and PE (31.6%). Remarkably, 15 out of 19 lipids (79%) contained at least one FA with 18 carbon atoms.

Lipids of clusters 4 and 5 (*n* = 198 and 312, respectively) were generally characterized by a strong concentration increase during adipogenesis. However, cluster 4 species increased substantially only until days 12–16 and decreased thereafter to half their maximum concentrations, while cluster 5 species reached a plateau at day 16. In cluster 4, TAG was the most abundant lipid class (75.8%), followed by PE (13.1%), and DAG (6.1%), while CE and all ceramide classes (CER, HCER, and LCER) were completely absent. The species of this cluster had the lowest number of C-atoms and DB within the dataset (Figure 6B–C). In cluster 5, the relative amount of TAG was even higher (88%) than in cluster 4 (Figure 6A). Owing to the high number of TAG species, clusters 4 and 5 also showed the highest lipid concentration values among all clusters (Figure 2B). Therefore, clusters 4 and 5 can be considered as "TAG clusters."

The lipids in cluster 6 (*n* = 56) exhibited a fluctuating averaged concentration profile throughout the investigated time period. Lipid species of all classes except LCER were represented in this cluster, which was in general rather heterogeneous in terms of composition (Figure 6A).

The correlations between the clusters might reveal new insights into lipid remodeling. In general, we observed positive correlations between clusters 1 and 2 as well as between clusters 4 and 5. Remarkably, lipids from clusters 1 and 2 were strongly negatively correlated with lipids in clusters 4 and 5, which might be an indicator of lipid remodeling. Specifically, many TAG species from cluster 5 had strong negative correlations with more than a dozen PE species of cluster 1 (between −0.7 and −0.96). These PE species mostly contained polyunsaturated and VLCFA, whereas the TAG species were carrying at least one LCFA (Table S5).

Moreover, some sphingomyelins, especially SM 20:0 (cluster 4), SM 22:0 (cluster 5), and SM 24:0 (cluster 5), had high negative Spearman's correlation coefficients (down to −0.84) with species from classes CER, HCER, and LCER (all in clusters 1 and 2, Table S5). These results point to regulatory interactions between the lipid species over several lipid classes.

In addition, we also found high positive correlation coefficients (mostly between 0.85 and 0.98) between DAG and TAG species. Those DAG and TAG with very high co-correlations mostly contained one or two of the most abundant fatty acids (C16:0, 16:1, 18:0, and 18:1) as side chains. Additionally, these TAG were characterized by total C-atom numbers between 42 and 54, which further implies that LCFA were the main constituents.

**Figure 6.** Clusters differed widely in their lipid class and lipid species composition. Panel (**A**) highlights the cluster compositions by showing the relative numerical proportions of the lipid class members. Panels (**B**–**D**) show the analysis of the FA side chain lengths and the number of DB based on the number of side chains of the lipid classes. Panel (**B**) shows the analysis of the class with three bound FA side chains, panel C that of lipid classes with two bound FA side chains, and panel D that of classes with only one FA side chain. Lipid species with patterns of decreasing concentration throughout adipogenesis had in general the highest numbers of DB and chain lengths, independent of their number of side chains.

#### **3. Discussion**

In the present study we characterized the different stages of human adipogenesis by using the Lipidyzer™ method [26] which enabled the simultaneous quantification of 743 lipid species of 11 different lipid classes in differentiating human SGBS cells. This global lipid analysis enabled us to identify correlations between lipid species over several classes and opened the possibility to generate hypotheses on lipid remodeling during adipogenesis.

Prior to the adipogenesis characterization, we analyzed the method performance in terms of linearity, repeatability, and background signals to be able to apply this novel methodology for accurate lipid analysis to SGBS cell culture samples. We used representative samples for undifferentiated and differentiated SGBS cells and validated the Lipidyzer™ method according to recently published studies that already showed good performance of the Lipidyzer™ method [21–25].

The lipid classes CE, CER, DAG, HCER, LCER, LPC, LPE, PC, PE, SM, and TAG could be analyzed with good linearities as indicated by coefficients of determination above 0.9 in at least one SGBS sample type. Furthermore, these lipid classes could be measured with high precision (CVs < 12%) which was determined by repeated analyses of QC samples from pooled SGBS homogenates. However, DCER and FFA had to be discarded from the data set, because the observed concentration levels were mostly in the range of the blank samples and therefore led to poor linearity and low precision. The detected high background signals of FFA in blank samples indicate high contaminations of these species that can be often found in glassware, pipette tips, and even in high-grade organic solvents [27]. We conclude that the Lipidyzer™ method was suited to reliably quantify 743 lipids from 11 lipid classes for the analysis of SGBS cell samples.

We are the first to use this novel method for the determination of lipid levels in differentiating human adipocytes. Additionally, we are also the first who characterized the human adipogenesis with this broad coverage of different lipid classes using only one method. Liaw and coworkers investigated the endpoints of adipogenesis (i.e., pre-adipocytes vs. fully differentiated adipocytes, day 12) in murine adipocytes of a large amount of lipids using an LC-MS/MSALL shotgun lipidomics approach [18]. They identified a shift from highly unsaturated VLCFA bound to the backbones of TAG, SM, cardiolipins, and ether-linked monoalkyldiacylglycerols in preadipocytes to more saturated LCFA in differentiated 3T3-L1 adipocytes. We could confirm these findings for VLCFA-carrying SM since we identified a decrease in their concentration levels until day 16. However, we measured an increase from day 16 to 20 almost reaching the concentration levels of day 0. The authors investigated adipogenesis at day 12, possibly missing the increase in the late phase of differentiation we observed. On the other hand, and different to the study of Liaw et al., we observed members of further lipid classes, namely CE, CER, DAG, HCER, LCER, LPC, PC, and PE species, having strongly decreased amounts of highly unsaturated VLCFA in fully differentiated adipocytes when compared to preadipocytes.

We followed the adipogenesis process over 20 days by analyzing samples from six time points (days 0, 4, 8, 12, 16, and 20). Successful cell differentiation of preadipocytes into lipid-laden adipocytes was confirmed on two different levels. First, we could observe the expected morphological changes during cell differentiation by microscopy. Second, the upregulation of *PPARG* and *CEBPA* expression showed strong transcriptional activation of cell differentiation. The very tight clustering of the data for samples from the same time point in the PCA and PLS-DA plots demonstrated generally high quality of the data obtained by the targeted lipidomics technique. Furthermore, the shift from differentiating (days 0–4) to maturating (days 4–12) SGBS cells became clearly visible by the change from PC1 to PC2 as the main contributor for cluster separation in the score plots. The observed shift coincides with findings from Halama and coworkers, who characterized murine adipogenesis of 3T3-L1 cells using a combined metabolomics and transcriptomics approach [16]. This shift during adipogenesis can be explained by the change from differentiation to maturation medium at day 4.

During ongoing adipogenesis we observed on the one hand strongly decreasing levels of CE, CER, HCER, and LCER, and on the other hand, substantially increasing levels of nearly all investigated glycerophospholipid classes, namely LPE, PC, PE, SM, as well as TAG. However, the temporal concentration courses of the individual lipid classes followed a particular pattern that was dependent on the developmental phases of the cells.

During the differentiation phase, we observed increasing concentrations of CER until day 4. The ceramides are known to be involved in signaling activity in cell cycle arrest and the inhibition of

cell proliferation during early adipogenesis [28]. Both processes are required for the induction of cell differentiation of preadipocytes into adipocytes [29]. Thus, the time courses of CER concentrations give rise to the hypothesis that at least some of these compounds were involved during cell differentiation signaling. After day 4, when cell differentiation was completed and maturation started, we observed decreasing levels of CER together with HCER and LCER accompanied by a simultaneous increase of SM. This could be explained by the remodeling of all ceramides into SM. This is further supported by the fact that we observed high negative Spearman's correlation coefficients between several SM species and CER, HCER, and LCER from the sphingolipid and PUFA cluster (cluster 2).

The differentiation phase was also characterized by decreasing concentration levels of CE. The CE operate as transport intermediates of cholesterol, which is an important component of the cell membrane [30,31]. Decreasing CE concentration levels indicate the release of cholesterol and its insertion into the cell membranes [32]. Incorporation of cholesterol decreases the flexibility of plasma membranes and thereby enables the morphological changes of the membranes that are important for differentiation and maturation. Furthermore, the cholesterol might also have been incorporated into triglyceride lipid droplet surfaces, serving as an intracellular free cholesterol reservoir [31]. These hypotheses are supported by negative correlations down to −0.78 between CE in the PUFA clusters 1 and 2 and TAG in the TAG clusters 4 and 5.

The second phase of adipogenesis, the maturation phase, was characterized by a strong increase of TAG from micromolar to millimolar levels. It is not surprising that TAG became the dominant lipid constituent of the adipocytes, as this reflects the function of adipocytes as sites for the storage and supply of fatty acids. The biosynthesis of TAG from DAG is corroborated by high positive correlations between these two lipid classes [33,34]. We observed a biphasic pattern of DAG and TAG concentration courses during adipogenesis. After a lag phase until day 4, considerable production of TAG via DAG reached a maximum at day 16 and then slowly declined. We conclude that the massive synthesis of TAG and their precursors started only after full differentiation from preadipocytes to adipocytes. As Collins and coworkers have shown by using 13C-labeled substrates, the massive generation of TAG is presumably based on de novo lipogenesis from glucose provided in the cell culture medium [19]. However, our correlation analysis revealed strong negative correlations between several PE species containing VLCFA and TAG. This might be an indicator of a possible contribution of PE containing VLCFA as an additional source for TAG synthesis during adipogenesis. The PE species might have been catabolized by phospholipase C to DAG, re-esterified to TAG, and incorporated into lipid droplets [35,36]. The increase of intracellular lipid depots is accompanied by expansion of the cell surface and volume, which requires larger amounts of the major membrane lipid classes like PC, PE, SM, and cholesterol [37–39]. Indeed, we observed simultaneous increases of PC, PE, and SM after day 4. Furthermore, LPC and LPE, which can be regarded as metabolic intermediates of PC and PE, also increased during adipogenesis [40].

Surprisingly, we observed the odd chain fatty acids C15:0 and C17:0 in CE, DAG, LPC, PC, PE, and TAG at non-negligible concentrations. Both fatty acids showed strongly increasing concentration time courses during adipogenesis. Roberts et al. also quantified increased odd chain fatty acid levels during cell differentiation of murine 3T3-L1 adipocytes [9]. We speculate that the increasing levels during adipogenesis might be explained by sequential peroxisomal fatty acid α-oxidation, which has been shown to occur in differentiating adipocytes [41].

Furthermore, some of the time courses of certain lipids might be explained by influences of the composition of the culture medium. To investigate a possible contribution of FA from the culture medium, we measured the FBS-containing medium for the cultivation of preadipocytes as well as the differentiation medium with the Lipidyzer™ method. Lipids with very long-chain PUFAs were highly concentrated in the FBS-supplemented medium compared with the levels in the differentiation medium lacking FBS (Figure S7). Therefore, we hypothesize that the measured concentration profiles of these lipid species during adipogenesis could have been artificially influenced by the cultivation with FBS-containing medium before the start of differentiation. As a result of this, the decreasing concentration levels of the VLCFA might also be explained by a lack of supply of these FA in the FBS-free differentiation and maturation media during adipogenesis.

It is important to keep in mind that cell culture experiments reflect artificial conditions. First, in vitro experiments often require high concentrations of growth factors, hormones, or several stimulation factors, which do not reflect physiological in vivo conditions. For instance, the differentiation of SGBS and other pre-adipocyte cells requires the corticosteroid dexamethasone and the insulin sensitizer rosiglitazone, which might strongly influence the lipidome. Indeed, Jeucken and Breuwens recently showed that rosiglitazone has an effect on the lipidome of HeLa cells [42]. In addition, the PPARγ agonist rosiglitazone as well as the endogenous myokine irisin can induce browning of white adipocytes [43]. Some recently published manuscripts showed this propensity of the SGBS cells towards a beige phenotype [43–46]. The used protocol for the SGBS cells requires an initial four-day stimulation with rosiglitazone for the induction of differentiation. This short time period might have an influence on the lipidome. However, undifferentiated SGBS cells behave very similar to human primary preadipocytes and the fully differentiated cells cannot be morphological distinguished from human primary adipocytes [47]. Moreover, one study compared SGBS cells, derived from subcutaneous adipose tissue of a male infant, with primary subcutaneous adipocytes from obese female patients [44]. The different confounders, obesity and sex, might also have a significant influence on the differentiation capacity and therefore also on the comparison of the two cell models. In addition, SGBS cells carry an FTO risk allele and the cells do not have a Simpson-Golabi-Behmel syndrome typical mutation in the glypican-3-gene (GPC3) what the name of the cell strain would suggest [47]. Nevertheless, future efforts will be necessary to confirm our findings in human primary subcutaneous adipocytes. Second, the lipid synthesis in adipose tissue *in vivo* is based not only on de novo synthesis from mostly glucose but also on circulating fatty acids in the bloodstream [48]. To sum up, cell culture experiments are indeed helpful to shed light on several cellular processes separately. However, owing to their simplicity and artificial nutrition, they cannot represent physiological in vivo conditions.

It also has to be mentioned that the Lipidyzer™ technology has some limitations. First, it was developed for the analysis of human plasma, which of course has a different lipid composition from (pre)adipocytes. The internal standard concentrations were therefore optimized for human plasma and may not match the actual situation in cell culture samples. Another issue concerns quantification. Although the Lipidyzer™ uses up to 10 internal standards (IS) per lipid class, there is no IS for each individual analyte. In addition, the method does not use a calibration curve for absolute quantification. Therefore, the measured absolute concentrations should be interpreted carefully. Furthermore, the Lipidyzer™ method is capable to determine the lipid species at the fatty acyl/alkyl level, except for the class of TAG. In case of TAG, the method cannot distinguish the *sn-1*, *sn-2*, and *sn-3* positions of the glycerol backbone and can as well not determine the exact positions of double bonds in the side chains. Besides, as the Lipidyzer™ method is a commercial assay with specialized and standardized software it is not possible to include other lipids into the method. This unfortunately limits the availability of further interesting lipid species such as for example signaling phospholipids.

Despite these limitations, we think that targeted analytical methods using multiple internal standards per lipid class—like the Lipidyzer™ technology—should be the methods of choice for the quantitative analysis of longitudinal samples with strongly different analyte concentrations and matrix conditions between sampling points. It has recently been shown by Chamberlain et al. that "due to the presence of matrix effects in untargeted, non-quantitative metabolomics, the signal intensity of any single analyte cannot be directly compared to the signal intensity of that same analyte (or any other analyte) between any two different matrices" [49]. This is of particular importance for ESI-MS-based lipid analytics, because matrix effects can vary considerably between lipid classes and even lipid species of the same class can respond differently to matrix effects depending on acyl chain length and degree of unsaturation [50,51]. Chamberlain et al. further concluded that "due to differences in ionization efficiency, the signal intensity of any single analyte cannot be directly compared to the signal intensity of any other analyte, even in the same matrix." Thus, any kind of correlation or network

analysis would be hampered with non-targeted or shotgun approaches. The application of IS can avoid or at least reduce the negative impact of matrix effects on the results. Non-targeted metabolomics or shotgun lipidomics approaches not including IS are more susceptible to matrix influences and should be therefore interpreted very carefully.

#### **4. Materials and Methods**

#### *4.1. Cell Culture, Cell Harvesting and Homogenization*

The Simpson Golabi Behmel syndrome (SGBS) preadipocyte cell strain was provided by Martin Wabitsch. The cells were cultivated and differentiated for 20 days, as described previously [52]. In brief, 50,000 preadipocytes per well were seeded in six-well plates in DMEM/F-12 medium (Thermo Fisher Scientific, Waltham, MA, USA), supplemented with 10% FBS (Biochrom, Berlin, Germany), 3.3 mM biotin (Merck, Darmstadt, Germany), and 1.7 mM pantothenate (Merck, Darmstadt, Germany) and grown at 37 ◦C and 5% CO2 in a humidified atmosphere. Cell differentiation was initiated when cells reached about 90% confluence. At that point, the medium was exchanged for serum-free medium supplemented with 10 μg/mL transferrin (Merck, Darmstadt, Germany), 0.2 nM triiodothyronine (T3; Merck, Darmstadt, Germany), 250 nM hydrocortisone (Merck, Darmstadt, Germany), 20 nM human insulin (Merck, Darmstadt, Germany), 25 nM dexamethasone (Merck, Darmstadt, Germany), 250 μM 3-isobutyl-1-methylxanthine (IBMX; Merck, Darmstadt, Germany), and 2 μM rosiglitazone (Biomol, Hamburg, Germany). After 4 days and then every fourth day, thereafter, the medium was replaced with serum-free medium containing 10 μg/mL transferrin, 0.2 nM T3, 250 nM hydrocortisone, and 20 nM human insulin (maturation medium). Cell morphology was monitored by microscopy.

Cell samples for quantitative real-time PCR (qRT-PCR) analyses were taken in quadruplicates (biological replicates) at each of 5 time points beginning with the start of differentiation (representing day 0), then on day 2, 4, 8, and 12 of adipogenesis. Cell samples were processed as described below. Cell samples for Lipidyzer™ analyses were taken at six time points beginning with the start of differentiation (representing day 0) and then on every fourth day of differentiation until day 20. At each time point, six samples (biological replicates) were collected. Harvesting, homogenization of cells, and normalization of measured metabolite concentrations to the cell number were performed as recently reported [53], with minor modifications. In brief, after one washing step with 6 mL of warm PBS per well of a six-well plate, the cells were scraped off the wells in 500 μL of extraction solvent of ice-cold 80% methanol per well using rubber-tipped cell scrapers (Sarstedt, Nümbrecht, Germany). Harvested cell-solvent suspensions of four wells were pooled into pre-cooled 2 mL microtubes (Sarstedt, Nümbrecht, Germany) containing 400 mg of glass beads (Bertin, Frankfurt, Germany). The samples were stored at −80 ◦C until further use. Homogenization of cells was performed immediately before analyses at 4–10 ◦C twice for 25 s at 5500 rpm using a Precellys24 (PeqLab, Erlangen, Germany). The resulting homogenates were used for lipidomics measurement by FIA-(DMS)-MS/MS as well as for DNA quantification (DNA content indirectly reflects the cell number of the sample and was determined for normalization purposes) [53].

#### *4.2. RNA Isolation and Quantitative Real-Time PCR (qRT-PCR)*

Total RNA of four independent biological replicates per group was extracted from cells using miRNeasy mini kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol. Synthesis of cDNA was performed by using the RevertAid First Strand cDNA Synthesis Kit (ThermoFisher Scientific, Dreieich, Germany) according to the manufacturer's instructions. Total RNA was reverse transcribed using an anchored oligo(dT)18 primer (5'-TTTTTTTTTTTTTTTTTTVN-3') in a final concentration of 0.5 μM for priming cDNA synthesis. For qRT-PCR, primers were designed using Primer3 to span at least one exon-intron boundary to avoid falsified amplification results [54]. Primers were synthesized by Metabion (Planegg, Germany) and sequences were as follows: PPARG\_for (5'-GACCACTCCCACTCCTTTGA-3'), PPARG\_rev (5'-GAGATGCAGGCTCCACTTTG-3'), CEBPA\_for

(5'-AACAGCTGAGCCGCGAACTG-3'), CEBPA\_rev (5'-CGGAATCTCCTAGTCCTGGCT-3'), TBP\_for (5'-CAGCCGTTCAGCAGTCAA-3'), TBP\_rev (5'-CTGCGGTACAATCCCAGAAC-3'). The amplification was performed on a QuantStudio Flex 7 Real-time PCR system (ThermoFisher Scientific, Dreieich, Germany) in triplicates using Power SYBR Green PCR Mastermix with ROX as passive reference (ThermoFisher Scientific, Dreieich, Germany) as follows: Denaturation at 95 ◦C for 10 min, 39 amplification and quantification cycles with 95 ◦C for 15 s and 60 ◦C for 1 min, and finally a melting curve program (95 ◦C for 15 s, followed by 60–95 ◦C with a heating rate of 0.1 ◦C/s) and continuous fluorescence measurement. The cycle threshold (CT) values were determined using the QuantStudio Flex 7 Real-time PCR system software. Relative gene expression was calculated using the comparative 2−ΔΔCT method [55]. Amplification efficiencies were determined based on the slope of the calibration curve consisting of five different cDNA concentrations each measured in triplicates and efficiencies were as follows: *PPARG* (106.3%), *CEBPA* (96.3%), and *TBP* (88.7%). The fold-change values for gene expression were normalized by the respective efficiencies using a published procedure [55]. Relative gene expression data for *PPARG* and *CEBPA* were subsequently normalized to the reference gene of tata-box binding protein (TBP; in pre-experiments tested to be suited) and the gene expression of samples at day 0 of adipogenesis.

#### *4.3. Hoechst Assay for DNA Quantification*

For DNA quantification, the fluorochrome Hoechst 33342 (ThermoFisher Scientific, Waltham, MA, USA) was diluted in PBS to a final concentration of 20 μg/mL. A total of 80 μL of this solution was pipetted into each of the wells of a black 96-well plate (F96; Nunc, Thermo Fisher, Schwerte, Germany). Then, 20 μL of vortexed cell homogenates or plain solvent (80% MeOH; blanks) was added to the Hoechst solution and mixed by pipetting. The plate was incubated in the dark for 30 min at room temperature. Fluorescence signals were read using a GloMax Multi Detection System (Promega, Mannheim, Germany), equipped with a UV filter (λex. = 365 nm; λem. = 410–460 nm, Promega, Mannheim, Germany) [53].

#### *4.4. Lipid Extraction and Targeted Lipidomics Analysis*

The Lipidyzer™ method (SCIEX, Darmstadt, Germany) was used to analyze the cellular lipidome. It detects lipids with fatty acid side chains of medium-chain (MCFA; C12), long-chain (LCFA; C13–C21), and very long-chain (VLCFA; C22–C26) lengths from 13 classes of lipids including cholesterol esters (CE), ceramides (CER), dihydroceramides (DCER), diacylglycerols (DAG), free fatty acids (FFA), hexosylceramides (HCER), lactosylceramides (LCER), lysophosphatidylcholines (LPC), lysophosphatidylethanolamines (LPE), phosphatidylcholines (PC), phosphatidylethanolamines (PE), sphingomyelins (SM), and triacylglycerols (TAG). The method allows the identification of lipid species at the fatty acyl/alkyl level (exception: TAG) [56]. The Lipidyzer™ method determines the sum of the numbers of C-atoms and double bonds (DB) for one fatty acid side chain as well as the sum of the C-atoms and DB of all three side chains. The notation rules from Liebisch and coworkers only know the case that either no fatty acid is known (e.g., TAG 52:2) or all three (e.g., TAG 16:0\_18:1\_18:1) [56]. Therefore, the nomenclature for TAG species in our study was adopted to these recommendations. The internal standard (IS) mixture (Avanti Polar Lipids, Inc., AL, USA) was prepared in accordance to the Lipidyzer™ manual. For QC samples, 250 μL of pooled cell homogenates were used, consisting in equal parts of undifferentiated, differentiating (day 8 of differentiation), and maturely differentiated cells (day 16). Three reference plasma samples (SCIEX, Darmstadt, Germany) of 100 μL in volume were spiked each with 50 μL of the QC spike mixture (SCIEX, Darmstadt, Germany) to investigate inter-run and inter-project effects. Lipids were extracted by two-phase separation using methyl *tert*-butyl ether (MTBE), methanol, and water [57]. Briefly, 250 μL of cell homogenates for the main experiments or 10–300 μL for method evaluation experiments, QC samples, or QC spiked plasma samples were transferred to 1.5 mL safe-lock reaction tubes (Eppendorf, Hamburg, Germany). For each time point, we took in total six biological independent cell samples. Next, 160 μL of MeOH

and 900 μL of MTBE were added to each tube and incubated for 30 min at 900 rpm and room temperature in a shaker. For phase separation, 500 μL of H2O was added to each tube, the mixtures were vortexed, and the tubes were centrifuged at 15,000× *g* for 4 min at RT. The upper organic phases were transferred into glass vials. The extraction step was repeated once and organic phases were combined. Organic solvents were evaporated to complete dryness under a stream of gaseous nitrogen and residuals were reconstituted with 250 μL of sample running buffer (10 mM ammonium acetate in dichloromethane:methanol (50:50 *v*/*v*)). Samples were then analyzed with the Lipidyzer™ method, consisting of a Sciex 5500 MS/MS QTRAP system (SCIEX, Darmstadt, Germany) equipped with a SelexION ion source for differential mobility spectrometry (DMS), in accordance with the manufacturer's instructions [24]. A sample volume of 50 μL was injected with a Shimadzu Nexera X2 liquid chromatography system (SCIEX, Darmstadt, Germany) at an isocratic flow rate of 7 μL/min. Data were acquired automatically with the Lipidyzer™ Workflow Manager software (version 1.0, SCIEX, Darmstadt, Germany). The obtained concentration values in nmol/g supplied by the software were converted to μmol/L with the assumption that 1 mL of cell culture sample was equal to 1 mL of plasma which is equal to1g[58]. Converted concentration values were normalized by Hoechst assay results.

#### *4.5. Data Analysis*

To trace the process of adipogenesis, the concentrations of the single lipid species, summed concentrations of lipid classes, concentrations of fatty acids, and percentage compositions were analyzed. Lipid species were completely excluded from the data set if concentration values were missing (NA) in more than 33.3% of the samples within a time point. Missing values were replaced by the respective minimal lipid species concentration measured, divided by <sup>√</sup>2 and multiplied by a randomly chosen factor between 0.75 and 1.25.

Statistical analyses and graphical illustrations of the lipidomic data were performed using the software MetaboAnalyst 4.0 [59], GraphPad Prism 8.1.1, and R 3.5.1 [60].

Univariate statistical analyses were performed using the Mann–Whitney U test and Kruskal–Wallis test with Dunn's post-hoc test. Spearman's rank correlation analysis was used to test correlations between lipid species. Prior to PCA and other multivariate statistical analyses, lipid concentrations were log-normalized and auto-scaled (mean-centered and divided by the standard deviation of each variable) to achieve a normal distribution of the data set. Averaged concentrations are shown with standard deviations.

#### **5. Conclusions**

We used the Lipidyzer™ technology to study the cellular lipidome during the development of preadipocytes into maturating and finally mature adipocytes in a human cell culture model. The switch from differentiating preadipocytes to maturating adipocytes became clearly visible at the lipidome level. The differentiation process was accompanied by increased concentrations of ceramides that are known to be involved in cell differentiation signaling. While these ceramide species decreased after completion of differentiation around day 4, massive lipid remodeling occurred during maturation of the adipocytes. This maturating phase was characterized by the substantial synthesis of DAG and TAG species. We furthermore observed increases of membrane lipids like PC, PE, and SM as well as their biosynthetic precursors. Moreover, we could also show that the compositions of the lipid species itself became more homogeneous during differentiation to highly concentrated saturated/monounsaturated LCFA with the four most abundant fatty acids being C16:0, C16:1, C18:0, and C18:1. Interestingly, VLCFA constantly decreased in almost all investigated lipid classes during the maturation process. High negative correlation coefficients between membrane lipids containing VLCFA and TAG species imply that these lipids might have served as additional sources for TAG synthesis. However, further studies are necessary to shed light on other lipid classes, lipid synthesis, degradation, and remodeling pathways during adipogenesis. For instance, fluxome-based approaches could be helpful to follow

especially polyunsaturated and VLCFA within the preadipocytes. Moreover, a multi-omics approach might be helpful to detect connections between different pathways within the lipidome as well as the connection with pathways of more polar metabolites. In addition, it would also be interesting to compare the lipidomes of SGBS cells and freshly isolated human adipocytes since the metabolism in these cells might better reflect the original situation in humans. We could also show that the cultivation of cells with FBS-containing medium might influence the metabolism of the cells from several days up to weeks later. We recommend analyzing the medium used for cell culture lipidomics/metabolomics studies to prevent misinterpretation of the data.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2218-1989/10/6/217/s1, Figure S1: Analytical evaluation of the Lipidyzer™ method for undifferentiated and differentiated cells revealed strong linearity for most of the analyzed lipid classes; Figure S2: Representative microscopic images of SGBS cells showed successful differentiation based on morphological changes and a strong increase in number and size of lipid droplets; Figure S3: Transcripts of the main adipogenic transcription factors PPARγ (*PPARG*) and C/EBPα (*CEBPA*) were highly upregulated during adipogenesis; Figure S4: Principal component analysis (PCA) score plot showing very clear clustering of lipid species regarding the different time points of adipogenesis; Figure S5: Fatty acid concentrations and compositions changed markedly during adipogenesis in all 11 lipid classes; Figure S6: Clusters from the spearman's rank correlation analysis consisted of distinct lipid class compositions; Figure S7: Polyunsaturated and very long-chain fatty acids were highly concentrated in FBS-containing medium used for cultivation before differentiation start; Table S1: Linearity testing of method; Table S2: Repeatability and background signal testing; Table S3: Top 30 variable importance in projection (VIP) of PLS-DA for component 1 and 2, respectively; Table S4: Significances of Kruskal Wallis test; Table S5: Spearman's rank correlation coefficients.

**Author Contributions:** Conceptualization, J.A. and M.H.; methodology, F.M. and J.T.; software, F.M., A.C., and M.H.; validation, F.M. and J.L.; formal analysis, F.M. and M.H.; investigation, F.M.; resources, J.A.; data curation, F.M. and A.C.; writing—original draft preparation, F.M.; writing—review and editing, F.M., G.M., J.T., J.L., M.W., J.A. and M.H.; visualization, F.M.; supervision, M.H.; project administration, J.A., J.T., and G.M.; funding acquisition, J.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** Part of this study was supported by EIT Health (eithealth.eu). EIT Health is supported by EIT, a body of the European Union.

**Acknowledgments:** The authors are grateful to Gabriele Zieglmeier for her technical assistance with the cell cultivation. Furthermore, we thank Maximilian Schmidtke for his support in RT-qPCR experiments.

**Conflicts of Interest:** The authors declare that they have no conflicts of interest with the contents of this article.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
