**Lipidomic Profiling of the Epidermis in a Mouse Model of Dermatitis Reveals Sexual Dimorphism and Changes in Lipid Composition before the Onset of Clinical Disease**

#### **Jackeline Franco 1, Bartek Rajwa 2,\*, Christina R. Ferreira 3, John P. Sundberg <sup>4</sup> and Harm HogenEsch 1,5,\***


Received: 16 June 2020; Accepted: 18 July 2020; Published: 21 July 2020

**Abstract:** Atopic dermatitis (AD) is a multifactorial disease associated with alterations in lipid composition and organization in the epidermis. Multiple variants of AD exist with different outcomes in response to therapies. The evaluation of disease progression and response to treatment are observational assessments with poor inter-observer agreement highlighting the need for molecular markers. SHARPIN-deficient mice (*Sharpincpdm*) spontaneously develop chronic proliferative dermatitis with features similar to AD in humans. To study the changes in the epidermal lipid-content during disease progression, we tested 72 epidermis samples from three groups (5-, 7-, and 10-weeks old) of *cpdm* mice and their WT littermates. An agnostic mass-spectrometry strategy for biomarker discovery termed multiple-reaction monitoring (MRM)-profiling was used to detect and monitor 1,030 lipid ions present in the epidermis samples. In order to select the most relevant ions, we utilized a two-tiered filter/wrapper feature-selection strategy. Lipid categories were compressed, and an elastic-net classifier was used to rank and identify the most predictive lipid categories for sex, phenotype, and disease stages of *cpdm* mice. The model accurately classified the samples based on phospholipids, cholesteryl esters, acylcarnitines, and sphingolipids, demonstrating that disease progression cannot be defined by one single lipid or lipid category.

**Keywords:** lipidomics; atopic dermatitis; SHARPIN-deficient mice; flow-injection mass-spectrometry; predictive elastic net

#### **1. Introduction**

Atopic dermatitis (AD) is a multifactorial inflammatory skin disease that affects people and domestic animals worldwide [1]. Multiple variants (endotypes) of AD occur based on differences in the genetic background of patients, environment, immune activation pathways, and epidermal barrier status [1–3]. The classical AD presentation includes increased IgE serum levels, increased concentration of type 2 cytokines [4,5], and filaggrin (*FLG*) mutations that underlie skin barrier dysfunction [6–8]. However, variants of AD with normal levels of serum IgE and an increase of Th22 and Th17 cytokines instead of type 2 cytokines also exist [7,9]. In addition, *FLG* mutations occur in only 10% to 30% of

AD patients [10,11]. The less common variants of AD may require different therapeutic approaches as standard forms of therapy could result in unsatisfactory outcomes. Currently, clinical assessment of disease severity and diagnosis of AD relies on subjective observation of clinical signs, which change with the chronicity of the disease phase [6,7]. Several assessment indices are used to diagnose and score the disease, but these have poor inter-observer agreement highlighting the need for molecular disease biomarkers [12–15].

Alterations in the skin lipid composition have been reported in AD patients regardless of the genetic background, immune response, and clinical presentation [16,17]. Investigation of the lipid composition of the stratum corneum of the skin across different analytical platforms revealed changes in ceramide (CER) structure and presence of shorter and more unsaturated free fatty acids (FFA) in AD patients compared to healthy subjects [18–21]. Others reported changes in the amounts of phospholipids (PL), cholesteryl esters (CE), and triacylglycerides (TAG) in atopic skin, sweat, and sebum compared with healthy controls [22–24]. Alterations in the lipid composition lead to a disorganized stratum corneum lipid matrix and impaired barrier function of the skin [18], which permits increased allergen penetration that induces or aggravates the inflammatory reaction [25,26]. The cause of these lipid changes is not well understood, and it remains uncertain whether they result from a primary defect or downregulation of lipid processing enzymes by type 2 cytokines released in the course of dermatitis [27,28].

*Sharpincpdm* mice (hereafter referred to as *cpdm* mice), which have a mutation that causes absence of the SHARPIN protein, develop a chronic proliferative dermatitis that is very similar to human AD. The condition is characterized by pruritus, alopecia, and thickening of the skin, as well as accumulation of eosinophils, mast cells, M2 macrophages, and increased expression of type 2 cytokines [29,30]. In a previous study, we identified specific changes in ceramides and fatty acids in the epidermis of female SHARPIN-deficient mice with chronic proliferative dermatitis using a novel accelerated mass spectrometry strategy, multiple reaction monitoring (MRM)-profiling [31]. As the severity of the dermatitis rapidly increases with age, *cpdm* mice present a suitable model to identify lipid changes in the skin before the onset of clinical signs of inflammation and during progression of the dermatitis.

Lipidomics allows the detection and identification of a large number of molecules in a high-throughput manner aimed at the identification of new biomarkers for diagnosis and disease progression as well as novel targets for treatment [32]. These systems biology approaches yield complicated, high-dimensional data that should not be analyzed using naive univariate statistical methods as they may produce a high false-positive rate when predicting and classifying phenotypes. Consequently, this data requires multivariate approaches [33,34].

Although predicting phenotype from lipidomic data can be performed using various machine learning approaches, the critical question asked by biologists searching for a mechanistic model is the meaning of the statistical prediction. The black box predictors may be entirely accurate, but they do not allow easy formation of post-classification hypotheses regarding the causal relationship between the employed features, and the produced prediction. On the other hand, ante-hoc explainable models such as regression-based approaches can be used not only for supervised classification but also for the identification of critically important covariates, which can be further studied in pursuit of a mechanistic model [35]. Therefore, feature selection and reduction employing methods such as elastic-net (ENET) regularized regression are beneficial for finding key predictive features in the rich biological data and for identifying potential biomarkers amid the vast number of responses produced by systems biology methodologies [36–39]. Here we report the postulated biomarkers of AD, delivered via a multi-tiered feature selection strategy that processed the data generated by MRM-profiling in order to characterize lipid changes in the skin before the onset of clinical signs, both at the level of lipid categories and individual lipids ions. The method was used to investigate the association of the identified features with disease progression in male and female *cpdm* mice and their age and sex-matched wild type (WT) littermates. The study identified alterations in lipid composition preceding the onset of clinical dermatitis and a subset of lipid ions predictive of the disease stage of each sample. Additionally, the data demonstrated that the epidermis of female and male mice had distinct lipid profiles and differed in the lipid changes associated with disease progression.

#### **2. Results**

#### *2.1. Association of Sex and Genotype to the Lipid Composition of the Mice Skin*

Epidermal samples (*n* = 72, 36 *cpdm,* and 36 WT) were monitored for the presence of 1030 lipid ions belonging to multiple lipid categories. First, the collected data were pre-processed as described in the Methods section by executing log-ratio transformations, followed by single decomposition value (SVD)-driven principal component analysis. The result was visualized in the compositional principal component (CPC) space.

The CPC projection clearly differentiated samples by sex with the first component explaining 30.2% of the data variance, whereas the second component accounting for 22.4% of the variance was mostly associated with the genotype (Figure 1). The list of transitions driving the separation of samples in the CPC score plot is provided in Table S1.

**Figure 1.** Monitored lipid ions in male and female *cpdm* and wild type (WT) epidermis by multiple reaction monitoring (MRM) scans in positive ion mode. Discrimination of the sex, as well as the genotypes of WT and *cpdm* mice (including non-lesional samples), was observed by compositional principal components (CPC) projection. (**A**) score plot of CPC analysis. (**B**) violin plots representing the separation of samples by sex and genotype. CPC 1 explained 30.2% of the variability of the data separating the samples by sex. CPC 2 explained 22.4% of the variance and was aligned with the genotype.

The visualization demonstrates that the information regarding the sex and genotype is encoded in the lipidomic profile of the sample. However, the CPC projection does not provide an actionable input from the perspective of feature selection or causal explanation. The 296 lipids in the top 25-percentile of accounted variance in CPC 1 contribute only between 0.167 and 0.31 percent to the representation. Similarly, for CPC 2, the individual contribution of each lipid in the top 25-percentile group ranges from 0.137 to 0.37. Therefore, a feature selection strategy is necessary.

The feature selection involved a two-tier selection, including a univariate step followed by the creation of a feature-ranking ENET regression model able to separate the samples into classes based on sex and genotype. The analysis was performed assuming a binary case (for sex and genotype data). These tasks were approached using two different methods: Analysis of CPC-compressed compositional features representing lipid categories, and analysis of individual ions regardless of the category.

#### 2.1.1. Selection of Predictive Lipid Categories for Sex and Genotype

The compressed features identified glycerolipids CPC 1, phospholipids CPC 1, and sphingolipids CPC 4 as most capable of separating samples by sex in the first selection step. For these features, the effect size expressed as η<sup>2</sup> ranged from 0.12 to 0.22 (Figure 2A). The η<sup>2</sup> of 0.22 is equivalent to Cohen's *f* = 0.53, which in a univariate model with two groups is equal to Cohen's *d* = 1.06, signifying a very substantial effect size. The subsequent feature-ranking ENET selected the sphingolipids CPC 4 and 5, phospholipids CPC 1, and glycerolipids CPC 1 as the most critical features to classify the samples by sex (Figure 2B). The classifier built using the 20-top composite CPC features had an overall accuracy of 0.76, CI0.95 = (0.74, 0.77).

The univariate selection of compressed features for the binary genotype classification (WT vs. *cpdm*) identified phospholipids CPC 3, glycerolipids CPC 5, cholesteryl esters CPC 2, acylcarnitine CPC 2, and sphingolipids CPC 3 and 1, as the most predictive. The observed η<sup>2</sup> ranged from 0.26 to 0.73 for the top features (Figure 2A). The subsequently trained ENET identified phospholipids CPC 3, glycerolipids CPC 5, sphingolipids CPC 3 and 1, and cholesteryl esters CPC 2 as the top features in terms of importance (Figure 2B) and the model approached 100% accuracy (CI0.95 from 1 to 0.95).

**Figure 2.** Importance of compressed categories for sex and genotype classification. (**A**) univariate linear model—selected lipid categories compositional principal component (CPC) based on their effect size (η2). (**B**) top 10 CPC of lipid categories ranked by the multivariate elastic net model for sex and genotype, which generate a prediction accuracy of 76% for the sex and 100% for the genotype.

These results demonstrate that the composition and abundance of the same lipid categories carry information regarding the sex and genotype status of the tested animals. Similar to the CPC visualization incorporating all the lipid ions simultaneously, we observed that the information regarding sex and genotype is present in multiple compositional principal components.

#### 2.1.2. Individual Lipid Ions Feature Selection for Genotype

Following the analysis of lipid categories, we performed a feature importance analysis for individual lipid ions. The study was conducted for genotypes (*cpdm* vs. WT) in a binary setting, filtering out the influence of sex. It is important to emphasize that even though this model used the individual lipid ion features, the tentative chemical attribution of the measured ions was not independently confirmed to eliminate the likelihood of isotopic interferences.

In the univariate step, we applied a linear model-based filter retaining only the features associated with genotype class (adjusted *p* < 0.01) and not strongly associated with sex (adjusted *p* > 0.05). The 100 lipid ions with the highest partial η<sup>2</sup> (ranging from 0.26 to 0.76) were selected for further analysis. In the top-100 group, acylcarnitines and phospholipids were by far the most prevalent. However, among the top ten features, there were nine phospholipids and one ion associated with sphingolipids. Table 1 shows the highest-scoring lipid ions as well as their weak effect size associated with sex.


**Table 1.** List of top 10 lipid ions ranked by the effect size of the univariate models linking genotype with the lipidomic profile.

\* Subject of possible isotopic interferences.

As in the previous analysis task, the second filtering step included an ENET regression used to filter and rank the lipids pre-selected by the univariate step. The trained ENET achieved an overall accuracy of 0.99, CI95% = (0.924, 1). The most predictive lipid ions are summarized in Table 2. Among the selected lipids, five were phospholipids, two glycerolipids, and one was identified as a sphingolipid.

**Table 2.** List of top lipid ions ranked by importance score for prediction of genotype using the elastic net model.


\* Subject of possible isotopic interferences.

#### *2.2. Selection of Features Associated with Disease Progression*

#### 2.2.1. Compositional Principal Component Analysis and Data Visualization

To study epidermal lipid changes associated with disease progression, a multiclass case was considered instead of a binary case. The *cpdm* mice were further divided into subclasses defined by the disease stage as non-lesional, established, and advanced. For general visualization of the data, we first computed CPC values using as input only the lipid data pre-selected in the previous binary step with the ENET filtering. The plot was prepared using the disease stage markings in a CPC space demonstrated that such a simple model was able to partially delineate the controls (independently of their age) and the levels of the *cpdm* genotype (Figure 3).

**Figure 3.** Lipid ions delineate disease stage groups. CPC analysis of all lipid ions plotted vs. disease progression. The model was able to delineate the controls and the three experimental groups of *cpdm* mice (η<sup>2</sup> = 0.86, *p*-value < 0.001).

#### 2.2.2. Compressed-Feature Selection for Disease Progression

The disease-progression analysis, performed in a univariate setting, pointed to phospholipids CPC 3, glycerolipids CPC 5, and cholesteryl esters CPC 2 as the most informative compressed features (Figure 4A). The top features associated with disease progression produced η<sup>2</sup> ranging from 0.63 to 0.8. It is important to note that the features predicting disease progression were the same as those that separated WT from the broad *cpdm* group containing animals in all disease stages. The feature selection and ranking task performed by the ENET again identified phospholipid CPC 3, cholesteryl esters CPC 2, and glycerolipids CPC 5, as the top features in terms of importance. Interestingly, the highly ranked features were not equally important for all the disease stages (Figure 4B).

The disease progression prediction with an ENET classifier using a multinomial model achieved an overall accuracy of over 0.81, CI95% = (0.71, 0.9). The substantial part of the observed inaccuracy was caused by the high similarity between samples from the adjacent "established" and "advanced" stages of the disease. This effect is also demonstrated by the difference between the unweighted and weighted Cohen's κ values (0.725 and 0.841, respectively).

**Figure 4.** Importance of compressed categories for disease progression classification. (**A**) univariate linear model—selected lipid categories CPC based on their effect size (n2). (**B**) contribution of the top 10 CPC to the classification of disease progression categories by the multivariate elastic net model with an accuracy of 81%.

2.2.3. Individual Lipid Ions Feature Selection for Disease Progression

The univariate feature selection step for disease progression selected ions with η<sup>2</sup> ranging from 0.22 to 0.73 and phospholipids dominated the very top of the list. The following multivariate analysis, performed by training an ENET, found a more diverse set of ions, some of them distinctly associated with a particular disease stage, but less useful for predicting others. It is an expected characteristic of a multivariate model, which combines all the features and their predictions to form a functional classifier. The ions contributing highly to the prediction of progression are listed in Table 3, and the results are illustrated in Figure 5.


**Table 3.** List of top 10 lipid ions ranked by importance score for prediction of disease progression in the elastic net model.

**Figure 5.** Epidermal lipid ions predictive of disease progression in mice. Representation of six lipids from the epidermis of WT and *cpdm* mice identified as predictive of disease stage in a sex-independent manner. Lipid features emphasize differences between controls and the various stages of the disease.

The ENET classifier trained on the disease progression data was able to classify the 36 *cpdm* and 36 WT samples into groups, including the control and the three disease stages with an overall accuracy of 0.79, CI95% = (0.67, 0.87) when classified using weighted classes and 0.95, CI95% = (0.88, 1) if the synthetic minority sampling technique (SMOTE) algorithm was used for correcting the class imbalance (Figure 6). The training used the variations transformed features corresponding to the presence of phosphatidylcholines, cholesteryl esters, acylcarnitines, and a glycerolipid-containing triacontanoic acid fatty acyl residue.

**Figure 6.** Classification of samples into disease progression groups. Parallel plot illustrating the classification of individual samples by elastic net regression. Each line represents a sample, and the highest point in the line corresponds to the group where the sample would be classified with higher probability.

#### **3. Discussion**

Lipids comprise a highly diverse group of molecules that play an essential role in the biology of the skin, and the relative proportions of different lipids are associated with the normal physiological functions of this organ [40–42]. In this study, we analyzed the relation between the lipids detected by MRM-profiling (lipidomic profile) and the observed genotypes using two machine learning feature-selection approaches. First, we computed a set of compressed features using compositional principal components to represent each of the lipid categories analyzed. These features easily separated male and female samples indicating a strong influence of sex on the epidermal lipid composition in mice. The first CPC summarizing variance in all the lipids was associated with clustering by sex, rather than by genotype. This result shows that the biochemical variability related to sex was dispersed among many lipids creating an effect more substantial than the one associated with the genotype. This result is in agreement with studies of skin-surface lipid clusters in humans where samples from males and females were distinguishable, but no significant difference between atopic or healthy subjects was observed [23,43]. However, it is necessary to note that the large variance visualized by CPC 1 (Figure 1) does not unequivocally demonstrate the importance of the differences between lipid composition in males and females, as it may also emerge from the fact that the sexual dichotomy has a high signal-to-noise ratio in the lipidomics data.

The multivariate analysis performed by the ENET demonstrated that classification to the male and female group was influenced mostly by sphingolipids (specifically, the compressed feature sets CPC 4 and CPC 5). The biological function of sphingolipids is determined by their composition, particularly the type of sphingoid base and the number of carbons and hydroxyl groups on the acyl chains, and their synthesis is affected by gonadal hormones in mice [44]. Several studies have shown alterations in ceramides, a sphingolipid, in the epidermis of AD patients [17,45–47]. However, conflicting results have been reported for changes in ceramides in non-lesional skin, probably because not all ceramide

species are altered at the same stage of the disease and/or by the same mechanisms in males and females [23,48]. Our results show that sexual dimorphism is related strongly to the relative amounts of epidermal lipids in mice, and suggest that sex-related differences in the lipid biology of AD should be further investigated as they may partially explain the contradictory results regarding changes in ceramides in AD patients' skin [47].

A comparison of lipid categories in the *cpdm* and WT phenotypes by either univariate method or ENET showed differences driven by phospholipids and glycerolipids. The alteration in the lipid composition of the epidermis is a hallmark of AD associated with impaired barrier function of the skin [41,49], but whether these changes are primary or caused by the inflammatory process remains elusive [50,51]. Lipidomic and transcriptomic analysis of atopic patients have shown a global alteration of fatty acids caused by the interrelationship of type 2 cytokines and lipid elongase enzymes [52]. Our analysis demonstrates that the presence of phospholipids and glycerolipids in the epidermis was altered before any clinical signs of disease in the skin of the *cpdm* mice. Phospholipids and glycerolipids are essential for cellular and subcellular membrane dynamics and share common metabolic intermediates [53]. Phospholipids can also be secreted in the lamellar bodies of the epidermis along with the enzymes that use them as a substrate for ceramide synthesis [54]. Alterations in their concentrations may affect skin barrier function, cell metabolism, and inflammatory cell signaling, as they can carry esterified fatty acids that generate lipid mediators of inflammation by undergoing fatty acyl remodeling [55]. Likewise, lipids resulting from an aberrant lipid metabolism may be incorporated into membranes as phospholipids are in constant flux [55,56]. In agreement with our results, changes in phospholipids were reported to be present in the serum of atopic patients compared with controls [57]. In another mouse model of AD, NC/Nga mice, phospholipids were decreased in plasma, and oral supplementation of plasmalogens increased the phospholipid concentration in the skin of the mice and improved the skin condition [58]. In human skin, one study reported an increase of phospholipid content in AD patients compared to healthy subjects [24], while others showed a global change in phospholipids with an increase in the presence of shorter acylated fatty acids [52]. In addition to species differences, the disparities in sample preparation approach (epidermis vs. total skin), variations in analytical methods, and identification of the individual lipid species may account for these inconsistent results.

Changes in different categories of lipids were associated with stages of the disease as the severity of dermatitis increased. The control epidermis could be discriminated from the *cpdm* samples mostly by phospholipids and glycerolipids. However, other lipid categories were required to separate the disease stages of the *cpdm* mice samples. Classification of stages was influenced by acylcarnitines, cholesteryl esters, and sphingolipids. Acylcarnitines are fundamental for β-oxidation by delivering fatty acids to mitochondria and peroxisomes as an energy source; therefore, their dysregulation could potentially redirect fatty acids towards increased biosynthesis of phospholipids and other lipids, rather than being used as a source of energy [59,60]. Acylcarnitines were found dysregulated along with the phospholipid content in the serum of atopic patients [57]. Acylcarnitines also play a role in inflammatory processes, and their accumulation has been linked to lipotoxicity that results in apoptosis [61]. Cholesteryl esters, like acylcarnitines, are carriers of fatty acids and serve to store cholesterol in lipid droplets for its transport [62,63]. Free cholesterol is an essential constituent of the epidermis [19]; however, cholesteryl esters are less polar than free cholesterol conferring the skin with enhanced hydrophobicity for the barrier function [64]. Increased levels of free cholesterol in the skin of atopic patients [51], along with reduced levels of cholesteryl esters associated with high-density lipoproteins [65] suggest alterations in systemic cholesterol homeostasis, as reported in cardiovascular diseases [66]. The data showing an influence of sphingolipids on the classification model reproduces our previous results in which ceramides were identified as necessary for the classification of samples from *cpdm* mice with advanced dermatitis [31]. Changes in ceramide content are correlated with impaired barrier function and increased transepidermal water loss [51,67]. They are also associated with the disease severity, resulting in considerable changes in lesional compared to non-lesional

skin [68]. Interestingly, sphingosine was an exception as its concentration was increased more in non-lesional than in lesional skin.

This study demonstrated an association between the presence of particular lipids in the epidermis and the occurrence of chronic proliferative dermatitis in mice. The link was displayed by a data-driven selection of MRM-extracted lipid features, training of a classifier, and identifying the key features contributing to the high classification accuracy. Although the presented procedure relies entirely on a statistical model, we hypothesize that the selected features will contribute to mechanistic insights if studied further. This argument is built upon the notion that the lipids providing high classification accuracy must be involved in the processes causing the phenotypic changes. Given that the lipid composition fingerprint differed not only between visibly lesional *cpdm* skin specimens and controls but was also altered in asymptomatic samples, the identified lipids could be considered candidates for predictive disease biomarkers. The current study involved an exploratory screening design that compares the lipid profiles of the groups using similar amounts of samples. The data analysis was performed employing compositional (relative) representation. Our study's scope did not include the validation of informative lipids by LC-MS/MS with the addition of internal standards; therefore, the isotopic and isobaric overlap may have occurred, and other lipids than the ones labeled with tentative attributions may have contributed to the reported predictive features. However, our previous work utilizing the relative amounts of ceramides demonstrated a concurrence between the relative values and the results obtained using quantitative LC-MS/MS [31]. Our analysis showed that not a single lipid (or lipid category) is modified sufficiently to be the sole differentiating factor. The accurate separation between healthy and diseased animals required an entire vector of lipid features, including phospholipids, acylcarnitines, cholesteryl esters, and sphingolipids. This result points to a multifaceted and multivariate nature of AD-associated lipid alterations in the skin.

The reliance on a classification tool to extract the most predictive lipids also defines the limitation of the presented approach. The task becomes particularly challenging when facing a severe class imbalance, as in the case of disease progression [69]. It is evident that sample availability determines the training performance, which in turn affects the robustness of the feature selection. We are aware of this limitation and hope that continuing research will allow for larger sample sizes, and correspondingly more confident analysis.

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

#### *4.1. Animals*

72 male and female C57BL/KaLawRij-*Sharpincpdm*/*Sharpincpdm* RijSunJ (*cpdm*) mice and WT littermates were obtained from the Jackson Laboratory and housed at 2 to 4 animals per box with food (Envigo) and water ad libitum. Room temperature was maintained at 20 ± 2 ◦C and relative humidity at 50 ± 15% with a 12/12-h light/dark cycle. Then, 18 WT males and 18 WT females and their *cpdm* littermates were divided into three groups of six males and six females with different ages and disease stages. The disease progression corresponded to non-lesional (5 weeks of age), established (7 weeks), and advanced (10 weeks) stages. Mice from the non-lesional group had no clinical signs of dermatitis on the dorsal or abdominal skin. The mice in the established group displayed erythema, moderate scaling, and mild alopecia of the dorsal and ventral skin. At 10 weeks, dermatitis covered most of the body with significant hair loss, erythema, thickening, and scaling. Mice were euthanized at 5, 7, or 10 weeks of age by CO2 asphyxiation and cervical dislocation. The animal experiments and procedures were conducted in accordance with the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. The protocol was approved by the Purdue University Animal Care and Use Committee (PACUC protocol 111001019).

#### *4.2. Epidermis Isolation and Lipid Extraction*

Sample collection and lipid extraction were performed as previously described [31]. Briefly, a 1 by 2 cm slice of dorsal skin was collected, and after incubation with Thermolysin (from *Geobacillus stearothermophilus*, Sigma-Aldrich, St. Louis, MO, USA) dissolved in HEPES buffer, the epidermis was peeled off and stored at −80 ◦C until extraction. Tissue was weighed and homogenized in 250 μL of ultra-pure water using Precellys24 tissue homogenizer (Bertin Technologies, Rockville, MD, USA). The homogenate was submitted to a Bligh and Dyer [70] liquid–liquid extraction, and the organic phase was collected and dried in a concentrator. Samples were resuspended in 40 μL of 3:1 (*v*/*v*) acetonitrile (ACN)/chloroform, then diluted 50× with ACN/methanol/ammonium acetate 300 mM at 3:6.65:0.35 volume ratio for mass spectrometry analysis.

#### *4.3. MRM-Profiling Method Development and Sample Screening*

A composite sample of each group from the testing set was created by pooling aliquots of 5 μL from each specimen in the group. The composite samples were analyzed using a previously described methodology of MRM-profiling discovery experiments [31]. Briefly, neutral loss (NL) and precursor ion (Prec) scans were used to profile phospholipids, acylcarnitines (AC), sulfatides, cholesteryl esters, ceramides, glycerolipids with diverse fatty acid acyl residues, triacylglycerides, and free fatty acids in positive and negative ion modes [31,71–75]. Using a micro-autosampler (G1367A), 8 μL of the sample was directly delivered into a QQQ6410 triple quadrupole mass spectrometer (Agilent Technologies, San Jose, CA, USA) equipped with an ESI ion source. A cap pump (G1376A) was used to flow acetonitrile plus 0.1% formic acid at a rate of 5 μL/min. The source capillary and multiplier voltages were 3500 V and 300 V, respectively. The collision energy voltage was 2 V for the negative ion mode methods. In positive ion mode, the collision energies varied according to the lipid classes. For ceramides, phosphatidylethanolamines (PE), and lipids with arachidonate acyl residue and oleate acyl residue, the collision energy was set at 22 V, for phosphatidylcholines and sphingomyelins (SM) at 20 V, for phosphatidylserines (PS) and phosphatidylinositols (PI) at 16 V, for CE at 17 V and for acylcarnitines the collision energy was set at 30 V. The fragmentation voltage of all the methods was 100 V. In total, 80 different discovery scans were performed, producing 1030 informative lipid ions. The parent and the fragment were collected and organized as transitions in 6 different methods of 2 min each (Table S2). The individual samples were flow-injected 6 times to cover all the monitored lipid ions. The raw data files are deposited in the public proteomics repository MassIVE (http://massive.ucsd.edu) using the identifier: MSV000083884. The tentative identification of lipid ions was performed through MS/MS experiments and by using reference databases, such as the Lipid Maps database (http://www.lipidmaps.org/) and METLIN (https://metlin.scripps.edu). Validation of the method by liquid chromatography-mass spectrometry has been previously reported along with the linearity and dynamic range of over four orders of magnitude from 1 to 10,000 ppm [31].

#### *4.4. Data Analysis*

Using MSConvert (http://proteowizard.sourceforge.net), the files were converted into the mzML open-source format, and an in-house script was used to obtain the ion intensities of each *m*/*z* monitored. The relative amounts of each *m*/*z* were used for data analysis.

The visualization and subsequent selection of lipid categories /groups and individual lipid ions associated with the *cpdm* genotype were performed following normalization of the signals to 1 using the total ion count and subsequent transformation using the isometric log-ratio function or centered log-ratio function. For visualization of all the overall characteristics of the data, the pre-processed input was compressed using singular value decomposition. The resultant first two compositional principal components were employed to illustrate the tendency of the samples to separate themselves in the reduced dimensionality space into clusters according to the sex, the genotype, and the disease severity [76].

#### 4.4.1. Selection of Predictive Lipid Categories

To identify the categories of lipids associated with the sex or genotype, the pre-processing and compression were followed by a two-tier selection including a univariate step, and a multivariate step driven by ENET regression.

In the beginning, the measured lipid ions were annotated and assigned to one of the following categories: (1) acylcarnitine, (2) acylcarnitine or glycerolipids, (3) cholesteryl esters, (4) DAG, (5) glycerolipids, (6) phospholipids, (7) phospholipids or cholesteryl esters, (8) phospholipids or glycerolipids, (9) sphingolipids, and (10) sphingolipids or glycerolipids. The overlap in categories reflects the uncertainty of the attribution due to the use of only one MRM related to a lipid candidate. Each of the sub-dataset was compressed using SVD to create compressed features sets. The number of retained columns (and by extension, the number of principal components used to represent the categories) was selected to retain 95% of the variance in each category. Therefore, the number of reduced composite features per class varied from 4 to 21. The resultant composite features for every lipid class were named "CPC", followed by the component number. For instance, "sphingolipids CPC 4" denotes the fourth principal component of the log-ratio transformed sphingolipids-class data.

It is important to emphasize that the SVD of the compositional data was not utilized here to enable a PCA-driven feature selection, but rather to produce a highly compressed input for the separate feature selection step. In other words, we are not claiming a direct association between the lipids that happen to display the most variation and the lipids (or lipid classes) that are most likely to be predictive and biologically significant.

The described data reduction process resulted in the creation of 57 compressed lipid-class features describing each of the samples. Subsequently, 57 linear models linking the computed features with sex, and another 57 models linking the features with the genotype (*cpdm* vs. control) were created. Finally, the third set of 57 linear models was computed to relate the features with the disease progression of *cpdm* mice (using the class assignment of control < non-lesional < established < advanced disease status). Benjamini–Hochberg *p*-value adjustment [77] was used to correct for false discovery. The features associated with genotype models having *p*-value < 0.05 (and η<sup>2</sup> effect sizes ranging from 0.73 to 0.1) were picked for the further feature selection step. For the sex-dependent changes in lipids, we also picked features with *p*-value < 0.05 (and η<sup>2</sup> effect sizes from 0.22 to 0.12).

#### 4.4.2. Feature Selection of Predictive Individual Lipid Ions

In a similar procedure, in order to recover the most predictive individual lipid ions (rather than lipid categories), we first created 1,030 linear models linking the log-ratio transformed relative amounts of every transition to genotype and sex. We pre-selected the features that might be associated with disease progression (either in sex-dependent or sex-independent manner) by selecting *p*-value < 0.01 for the criteria that were included and *p*-value > 0.05 for the factors that were ruled out. The lipids present in linear models connecting significantly with disease progression after Benjamini–Hochberg *p*-value adjustment (*p* < 0.01), but not being significant for sex (*p* > 0.05) were selected as sex-independent predictive lipids. About 50 ions were selected as possibly predictive and represented η<sup>2</sup> effect sizes ranging from 0.724 to 0.29.

#### 4.4.3. Predictive Elastic Net Regression

The final but critical step of these feature selection procedures involved the use of ENET regression [78]. ENET was employed either as a binary (for sex and *cpdm* vs. WT separation) or a multiclass classifier. This regression approach includes LASSO *L*<sup>1</sup> and ridge *L*<sup>2</sup> penalty terms leading to a predictive model operating in a reduced dimensionality of the data produced by the MRM-profiling:

$$\boldsymbol{\beta} = \underset{\boldsymbol{\beta}}{\text{argmin}} \Big( \left\| \mathbf{y} - \mathbf{X}\boldsymbol{\beta} \right\|^2 + \lambda \left( (1 - \alpha) \|\boldsymbol{\beta}\|^2 / 2 + \alpha \|\boldsymbol{\beta}\|\_1 \right) \Big) \tag{1}$$

In the ENET formula above, the input matrix **X** consists of all the pre-selected measured lipid ions (or pre-selected composite lipid categories), the output vector **y** describes the stages of the disease, and α, λ ≥ 0 are tuning parameters. The penalties included in the mathematical model are in the β<sup>1</sup> term which generates a sparse model by shrinking some regression coefficients to zero and, in the β<sup>2</sup> term which removes the limitation on the number of selected variables but encourages grouping effect, allowing similar features to be selected together. The individual lipid ions or the composite lipid categories with the larger absolute value of β are considered to be more predictive. The ENET simplifies to ridge regression when α = 1 and to the LASSO regression when α = 0.

The ENET regression was trained using the leave-one-out approach. Due to significant data imbalance, we used class weights (imposing a lesser penalty for errors in the majority class) or the SMOTE approach during training [79]. The resultant classifier allowed us to rank the lipid ions in terms of importance (ability to influence ENET prediction) using the absolute value of the non-zero coefficients.

To visualize the changes in the selected features, a central log-ratio transformation followed by standardization to the female WT subgroup was performed. Therefore, the y-axis in the figures shows the difference in relative lipid abundance as the number of standard deviations away from the WT-female group. The multiclass ENET prediction was illustrated using a parallel plot.

The statistical analyses were performed using R-language for statistical computing.

#### **5. Conclusions**

In this study, we paired an exploratory high-throughput lipidomics technique with rigorous machine learning analysis to rapidly screen for potential biomarkers in a mouse model of dermatitis. The measurements were performed using flow injection to the ion source of a triple quadrupole mass spectrometer, providing highly sensitive, but low-resolution mass data. The exploratory approach relied on product ions and neutral losses expected to be specific to the lipid classes, but not individual lipids; therefore, the detected lipids are assigned only tentative attributions. The approach revealed sexual dimorphism in the epidermal lipid profile, which was distributed throughout the lipid categories and identified sphingolipids as the best predictors for sex classification. Furthermore, epidermal lipid analysis allowed accurate classification of samples, not only by the genotype of the mice, cpdm vs. WT, but by the stages of disease progression. A panel of lipids comprised of phospholipids, acylcarnitines, sphingolipids, and cholesteryl esters was necessary to achieve successful classification into the different disease stages, showing that a single lipid or lipid category was not altered sufficiently to be the sole classifier. These results highlight the need to consider sex-related differences in the pathobiology of AD and the importance of building lipid panels that include lipids from different categories when investigating predictive biomarkers for AD.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2218-1989/10/7/299/s1, Table S1: Transitions for compositional principal component (CPC) score plot, Table S2: Compiled method for individual analysis of samples.

**Author Contributions:** Conceptualization, J.F., B.R. and H.H.; data curation, J.F. and B.R.; formal analysis, J.F. and B.R.; funding acquisition, J.F., J.P.S., and H.H.; investigation, J.F. and C.R.F.; methodology, J.F., B.R., C.R.F.; validation, J.F. and B.R.; project administration, H.H.; resources, J.P.S.; supervision, B.R. and H.H.; visualization, J.F.; writing (original draft), J.F., B.R. and H.H.; writing (review and editing), J.F., B.R., C.R.F., J.P.S. and H.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported in part by grants from the National Institutes of Health (AR049288) and the Purdue Institute of Inflammation, Immunology, and Infectious Disease (PI4D). JF was supported in part by a Colciencias fellowship.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **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* **LC–MS Lipidomics: Exploiting a Simple High-Throughput Method for the Comprehensive Extraction of Lipids in a Ruminant Fat Dose-Response Study**

#### **Benjamin Jenkins 1, Martin Ronis <sup>2</sup> and Albert Koulman 1,\***


Received: 22 June 2020; Accepted: 16 July 2020; Published: 17 July 2020

**Abstract:** Typical lipidomics methods incorporate a liquid–liquid extraction with LC–MS quantitation; however, the classic sample extraction methods are not high-throughput and do not perform well at extracting the full range of lipids especially, the relatively polar species (e.g., acyl-carnitines and glycosphingolipids). In this manuscript, we present a novel sample extraction protocol, which produces a single phase supernatant suitable for high-throughput applications that offers greater performance in extracting lipids across the full spectrum of species. We applied this lipidomics pipeline to a ruminant fat dose–response study to initially compare and validate the different extraction protocols but also to investigate complex lipid biomarkers of ruminant fat intake (adjoining onto simple odd chain fatty acid correlations). We have found 100 lipids species with a strong correlation with ruminant fat intake. This novel sample extraction along with the LC–MS pipeline have shown to be sensitive, robust and hugely informative (>450 lipids species semi-quantified): with a sample preparation throughput of over 100 tissue samples per day and an estimated ~1000 biological fluid samples per day. Thus, this work facilitating both the epidemiological involvement of ruminant fat, research into odd chain lipids and also streamlining the field of lipidomics (both by sample preparation methods and data presentation).

**Keywords:** odd chain lipids; lipid profiling; Folch; protein precipitation; sample preparation; relative lipid composition (Mol%)

#### **1. Introduction**

Lipids are generally understood as a class of molecules that have a high solubility in organic solvents and typically contain or originate from fatty acids. Although, lipids may be commonly derived, research has shown that there is a huge variety both structurally and functionally (potentially >40,000 [1]); where they play a vital role in energy production and storage [2,3], regulation and signalling [4,5], provide structure and support and membrane formation [6]. Lipids are now emerging as biomarkers of dietary/nutritional intakes [7] as well as indicators of pathophysiological status [8–11]. As a reaction of lipid-pathophysiological involvements, the field of lipidomics has emerged as a discipline that examines and quantifies a large proportion of the lipids present in a given sample set.

Lipidomics requires an effective isolation protocol that comprehensively extracts lipids from the sample as well as an analytical method that allows their identification and quantitation. The typical analyte isolation protocols (with/without minor adaptations) that are often used in the literature include three different liquid–liquid extractions: Folch and colleagues [12] (cited >65,000 times), Bligh and Dyer [13] (cited >52,000 times) or Matyash and colleagues [14] (cited >1000 times). Although these extraction protocols are heavily cited and do result in adequate results, there are several caveats with their use. Firstly, there is the need to perform duplicate extraction in non-fluid samples to ensure optimal recovery of the lipid analytes. This is extremely time consuming, especially for the Folch and the Bligh and Dyer methods. Secondly, there are reasonable concerns that using a biphasic extraction (producing immiscible aqueous and organic phases) may result in a loss of relatively polar lipids (e.g., acyl-carnitines and gangliosides) into the disposed aqueous fractions (consisting of mostly methanol and water in these extraction protocols: Folch and the Bligh and Dyer). There are publications that use a single phase extraction protocol but they do not appear to solve the problem of extracting the relatively more polar lipids since a mixture of methanol, chloroform and tert-butyl methyl ether were used [15,16].

The technique overwhelmingly used tor analysing the lipidome is mass spectrometry hyphenated with chromatography (LC–MS) due to its sensitivity and selectivity; furthermore by using a high-resolution accurate mass instrument (e.g., Orbitrap or Time-of-Flight instruments), a huge number of analytes can be analysed simultaneously. Reversed phase chromatography is the predominant chromatographic technique employed to separate the analytes before entering the mass spectrometer to determine their structure and concentration. Variants of a liquid chromatography method utilising a C18-column with a water and acetonitrile mix for the weak eluting mobile phase and acetonitrile and propan-2-ol for the strong eluting mobile phase are the most commonly used [17–22]. These reversed phase C18-column methods both separate lipid based on their lipid class assignment (i.e., either phosphatidylcholines or phosphatidyethanolamines head group) and their fatty-acyl composition (i.e., chain length and degree of unsaturation) with some degree of isomeric separation.

In this study, we present a lipidomics pipeline that including a novel analyte isolation protocol utilising a single phase, which results in a comprehensive lipid extraction suitable for a full range of lipid polarities (from polar to non-polar lipids species). This lipidomics method was tested, validated and then applied in a rat model investigating ruminant fat biomarkers via a beef tallow dose response dietary investigation.

#### **2. Results**

This lipidomics LC–MS method incorporating both of the described sample preparation protocols: protein precipitation (chloroform: methanol: acetone, ~7:3:4) and Folch liquid–liquid (chloroform: methanol: water, ~7:3:4), were utilised for the quantitation of lipids in liver samples from Sprague–Dawley rats who received one of four experimental diets overfed at 17% above matched growth.

A comparison of the two sample preparation methods on the extraction of the stable isotope-labelled internal standards are shown in the figure below (see Figure 1). A comparison on the samples' endogenous individual lipid classes are shown in the Supplementary Materials (see Supplementary Figure S1 and Table S1).

**Figure 1.** This figure shows the comparison between the two lipid extraction techniques regarding their extraction efficiency of the stable isotope internal standards from the rat liver samples (Folch liquid–liquid extraction: chloroform: methanol: water, ~7:3:4 , and Protein precipitation liquid extraction: chloroform: methanol: acetone, ~7:3:4 ). n = 34 rat liver samples per extraction method. The intensity of the internal standards were measured by liquid chromatography with mass spectrometry. The significance of the difference between the two extraction protocols are shown by the p-value star system; where *p* ≤ 0.05 was considered statistically significant (\* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001). Error bars represent ± standard deviation. Lipid internal standard include: Butyryl-d7-L-carnitine (abbreviated to IS\_Car\_4:0-d7), N-tetradecylphosphocholine-d42 (abbreviated to IS\_LPC\_14:0-d42), hexadecanoyl-L-carnitine-d3 (abbreviated to IS\_Car\_16:0-d3), heptadecanoic-d33 acid (abbreviated to IS\_FA\_17:0-d33), 1,2-dimyristoyl-d54-sn-glycero-3-[phospho-L-serine] (abbreviated to IS\_PS\_28:0-d54), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphoinositol (abbreviated to IS\_PI\_34:1-d31), N-palmitoyl-d31-D-erythro-sphingosylphosphorylcholine (abbreviated to IS\_SM\_34:1-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-[phospho-rac-(1-glycerol)] (abbreviated to IS\_PG\_34:1-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphate (abbreviated to IS\_PA\_34:1-d31), N-palmitoyl-d31-D-erythro-sphingosine (abbreviated to IS\_Cer\_16:0-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphocholine (abbreviated to IS\_PC\_34:1-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphoethanolamine (abbreviated to IS\_PE\_34:1-d31), glyceryl tri(pentadecanoate-d29) (abbreviated to IS\_TG\_45:0-d87).

As shown the intensity of seven of the internal standards are statistically significantly higher in the protein precipitation protocol when compared to the Folch liquid–liquid protocol (between ~30% to ~2500% higher), whereas, only two of the internal standards were higher in the Folch liquid–liquid protocol (between ~20% and ~37% higher). Additionally, there is far less variation in the protein precipitation liquid extraction protocol: 10 out of 13 internal standards had ~2% to ~72% less variation in their coefficient of variation (CV). The protein precipitation liquid extraction protocol

has also been shown to produce a significantly higher detection of the samples' endogenous lipids, both producing a higher total number of lipid detected (Folch-LLE: 455 lipid species, PPLE: 472 lipid species) and a statistically significantly higher total intensity for twelve of the sixteen lipid classes detected (see Supplementary Figure S1 and Table S1). Taken as a whole, the protein precipitation protocol (chloroform: methanol: acetone, ~7:3:4) showed a greater extraction capability across the full lipid hydro-philicity/phobicity range and across the internal standards. Additionally, the high throughput of the protein precipitation protocol allows for over 100 tissue samples to be extracted per day (including dissection, weighing and tissue extraction ready for LC–MS analysis), whereas, the Folch liquid–liquid protocol could take up to three- to four-times longer due to the necessity of duplicate extractions and the delicacy of liquid–liquid phase separation. The throughput of the protein precipitation protocol (chloroform: methanol: acetone, ~7:3:4) on fluid samples allows an estimated ~1000 biological fluid sample extractions per day (including aliquoting, sample extraction ready for LC–MS analysis) when utilising basic laboratory fluid handling equipment/robots (throughput data not shown here).

The liver lipid concentration (nM/mg) for each experimental diet group of rats are shown in the table below (see Table 1), along with the correlation (trendline equation, slope significance, R2 and successive change across the groups) of the measured lipid concentration with the percentage composition of ruminant fat (beef tallow) in each experimental diet. An R<sup>2</sup> threshold of 0.75, slope significance *p*-value < 0.05 and successive increase/decrease were set to establish if there was a strong correlation between the lipid concentration and the ruminant fat composition.








**Table 1.** *Cont*.














**Table** 







**TG\_(59:8)** TG\_(60:10)

**TG\_(60:12)** **TG\_(62:12)** **TG\_(62:13)**

 ± **9.09** ± **14.6**

 194 ± 194 **79.5** ± **74.4** **63.3** ± **43.7** **4.22** ± **2.92**

 73.7 ± 56.6

**39.6** ± **28** **31.3** ± **17.4** **2.86** ± **3.13**

 4.21 ± 4.15

**13** ± **12.1** **11** ± **5.05** **0.546** ± **1.15**

 0.838 ± 1.23 **1.42** ± **2.02** **0.863** ± **0.716** **0.0336** ± **0.0949**

**6.77** ± **11.9**

 ± **5.83** ± **9.99**

**1.75** ± **2.04**

 y = −0.0222x +

**y** = −**0.85x** + **12.4 \***

 y = −24x + 252

**y** = −**9.66x** + **107 \***

**y** = −**7.69x** + **85.4 \***

**y** = −**0.551x** + **6.13 \***

**Table 1.** *Cont*. **R2**

 0.89

 0.90

 0.12

 0.45

 1.00

 0.81

 1.00 **0.94**

 0.86 **0.94** **0.95** **0.95**

 Decreasing **Decreasing**

**Decreasing**

**Decreasing**

**Decreasing**

 Decreasing

 Decreasing

 Decreasing

**Successive Change**

**Across Groups**

#### **3. Discussion**

Out of 472 lipids detected and semi-quantified, 100 showed a strong relationship with the dietary intake of ruminant fat with 35 species increasing and 65 species decreasing as the percentage of ruminant fat in the diet increased (NB. ruminant fat as a percentage composition with corn-oil and medium chain triacylglyceride oil). Interestingly, ceramides generally increased, whilst cardiolipins, sphingomyelins and triacylglycerides generally decreased as the dietary composition of ruminant fat rose. According to the literature, a rise in liver ceramides is typically associated with aggravated non-alcoholic fatty liver disease (NAFLD) and insulin resistance [23], this in conjunction with a decrease in cardiolipins (which are indicative of mitochondrial remodelling and dysfunction [24]) may suggest that the changes in the experimental diets here are detrimental for these pathologies. However, there was a clear decrease in the triacylglycerides (particularly evident in the unsaturated odd chain triacylglycerides), which is explicitly representative of an ameliorated pathology [25]. As previously published [26], many NAFLD and insulin resistance factors were mitigated as the ruminant fat increased in these diets, including: a reduction in the total body weight (g), total fat mass (%), serum ALT (U/mL) and degree of steatosis determined by Oil Red O staining, notably, the inflammatory marker TNFα did not change significantly (trend: *p*-value = 0.52). A key characteristic of NAFLD development is the accumulation of hepatic triacylglycerides [27]; therefore, the data here suggests that these dietary changes may be beneficial for NAFLD and insulin resistance by aiding in a reduced hepatic triacylglyceride load: possible mechanisms here include a lower saturated fatty acid composition resulting in a lower fatty acid incorporation into hepatic triacylglycerides and/or a higher pass-through of the medium chain triglyceride oil directly into the mitochondria stimulation fatty acid metabolism [26]. Work presented by Gonzalez-Cantero and colleagues [28] showed that hepatic triacylglyceride content were correlated with insulin resistance and these relationships were independently to the inflammatory marker TNFα. Therefore, it appears that the hepatic triacylglyceride load may be paramount in the development of NAFLD and insulin resistance, which is supported in the literature [29].

According to the literature, odd chain fatty acids are considered biomarkers of their dietary intake and particularly accredited as a biomarkers of ruminant fat intake (e.g., milk, butter and beef tallow, etc.); however, there is a vast amount of conflicting data [2]. Some studies have shown both positive correlations (either individual odd chain lipids or total odd chain lipids) and some studies have shown there were no significant correlations. Although these studies may conflict in their findings they all present their data as relative compositions (Mol%), which is the typical way lipid data appear in the literature [30]. By expressing the lipid data as relative compositions (Mol%), it normalises the data to the total fat in that sample; however, presenting the lipid data in this way confounds the results by interconnecting the individual data points. This interconnection can cause false positive and/or false negative conclusions (type 1 and 2 errors), i.e., if a single lipid increases it will artificially decrease the other(s) due to the Mol% calculation. As shown in the figure below (see Figure 2), the concentration of the total lipids containing either even chain or odd chain fatty acids and a combination are shown. Although lipids containing odd chain fatty acids did increase, albeit not statistically significantly; *p*-value: 0.197, it was also not proportionate to the increase in dietary ruminant fat. As shown, the lipids containing even chain fatty acids did significantly inversely decrease (slope *p*-value: 0.0189) as the dietary ruminant fat increased. Interestingly, due to both the decrease in the even chain lipids and the consistency of the odd chain lipids, if the relative composition (Mol%) of the lipids were calculated, there appears to be a statistically significant increase in the odd chain lipids (see Figure 3); however, this is an artefact of changes in the even chain lipids and a consequence of interconnecting the data.

**Figure 2.** This figure shows the change in the liver lipid concentrations across the four high-fat diets fed to Sprague–Dawley rats (n = 8–9 per group): total odd chain lipids (symbol: **X**, trendline: , gradient: 0.183 <sup>±</sup> 0.0962, R<sup>2</sup> <sup>=</sup> 0.64, slope significance *<sup>p</sup>*-value: 0.197); total even chain lipids (symbol: , trendline: , gradient: <sup>−</sup>41.0 <sup>±</sup> 5.71, R2 <sup>=</sup> 0.963, slope significance *<sup>p</sup>*-value: 0.0189); total lipids containing both even and of odd chain (symbol: **<sup>O</sup>**, trendline: , gradient: <sup>−</sup>40.8 <sup>±</sup> 5.79, R<sup>2</sup> <sup>=</sup> 0.961, slope significance *p*-value: 0.0195. Lipid concentrations (μM/mg) are shown as means ± standard deviation and were extracted via the protein precipitation liquid extraction protocol (chloroform: methanol: acetone, ~7:3:4).

**Figure 3.** This figure shows the relative compositional (Mol%) change in the total odd chain lipids (symbol: **<sup>X</sup>**, trendline: , gradient: 0.0604 <sup>±</sup> 0.00203, R<sup>2</sup> <sup>=</sup> 0.998, slope significance *<sup>p</sup>*-value: 0.0011) and the total even chain lipids (symbol: , trendline: , gradient: <sup>−</sup>0.0604 <sup>±</sup> 0.00203, R2 <sup>=</sup> 0.998, slope significance *p*-value: 0.0011) across the four high-fat diets in Sprague–Dawley rats (n = 8–9 per group). Lipid compositions (Mol%) are shown as means ± standard deviation and were extracted via the protein precipitation liquid extraction protocol (chloroform: methanol: acetone, ~7:3:4). Diet one: 3.6% beef tallow; diet two: 6.3% beef tallow; diet three: 9.0% beef tallow; diet four: 11.7% beef tallow.

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

#### *4.1. Chemicals and Standards*

Stable isotope-labelled internal standards purchased from Sigma Aldrich (Haverhill, Suffolk, UK) include: N-palmitoyl-d31-D-erythro-sphingosine (abbreviated to IS\_Cer\_16:0-d31); order number: 868516P, 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphate (abbreviated to IS\_PA\_34:1-d31); order number: 860453P, 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphocholine (abbreviated to

IS\_PC\_34:1-d31); order number: 860399P, 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphoethanolamine (abbreviated to IS\_PE\_34:1-d31); order number: 860374P, 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-[phospho-rac-(1-glycerol)] (abbreviated to IS\_PG\_34:1-d31); order number: 860384P, 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphoinositol (abbreviated to IS\_PI\_34:1-d31); order number: 860042P, 1,2-dimyristoyl-d54-sn-glycero-3-[phospho-L-serine] (abbreviated to IS\_PS\_28:0-d54); order number: 860401P, N-palmitoyl-d31-D-erythro-sphingosylphosphorylcholine (abbreviated to IS\_SM\_34:1-d31); order number: 868584P. Stable isotope-labelled internal standards purchased from QMX Laboratories Ltd. (QMX Laboratories Ltd., Thaxted, Essex, UK) include: Heptadecanoic-d33 acid (abbreviated to IS\_FA\_17:0-d33); order number: D-5261, N-tetradecylphosphocholine-d42 (abbreviated to IS\_LPC\_14:0-d42); order number: D-5885, Glyceryl tri(pentadecanoate-d29) (abbreviated to IS\_TG\_45:0-d87); order number: D-5265, Butyryl-d7-L-carnitine (abbreviated to IS\_Car\_4:0-d7); order number: D-7761, Hexadecanoyl-L-carnitine-d3 (abbreviated to IS\_Car\_16:0-d3); order number: D-6646.

Quality control standards (LIPID-QC) purchased from Cayman Chemical Company (Cambridge Bioscience, Cambridge, UK) include: Lysophosphatidylcholines (egg); order number: 24331, Phosphotidylcholines (egg); order number: 24343, Lysophosphatidylethanolamines (egg); order number: 25844, Phosphatidylethanolamines (bovine); order number: 16878, Phosphotidlethanolamine (soy); order number: 25845, Lysophosphatidyinositols (porcine liver); order number: 26016, Phosphatidylserines (soy); order number: 25847, Ceramides mixture; order number: 22853, Ceramides (non-hydroxy); order number: 24833, Ceramides (hydroxy); order number: 24834, Sphingomyelins (from bovine spinal cord); order number: 22674, Sphingomylins (egg); order number: 24345, Phosphatidylglycerols (egg); order number: 25846, Phosphatidic acid (egg); order number: 24344, Sulfatides (bovine); order number: 24323, Purified mixed gangliosides (bovine); order number: 24856, TLC Neutral Glycosphingolipid Mixture (bovine and porcine); order number: 1505, 2-Palmitoyl Glycerol; order number: CAY17882, 1,2-Dipalmitoyl-sn-glycerol; order number: CAY10008648. Quality control standards purchased from Sigma Aldrich include: Soy PC (95%); order number: 441601G, C18(Plasm)-18:1-PC; order number: 852467C, Brain CPE; order number: 860066P, Liver PI; order number: 840042P, Brain lyso PS; order number: 850092P, Milk SM Sphingomyelin (Milk, Bovine); order number: 860063P, Galactocerebrosides from bovine brain; order number: C4905, Glucosylceramide (Soy); order number: 131304P, Triglyceride mix, C2–C10; order number: 17810-1amp-s, Fish oil from menhaden; order number: F8020, Anhydrous butter fat, Cardiolipin solution from bovine heart; order number: C1649, Brain PI(4)P; order number: 840045P.

Commercially available blank human serum was purchased from BioIVT (Royston, Hertfordshire, UK; order number: HUMANSRMPNN. All solvents and additives were of HPLC grade or higher and purchased from Sigma Aldrich unless otherwise stated.

LIPID-IS: the lipid stable isotope-labelled internal standard was prepared by dissolving each of the individual lipid standards into chloroform: methanol (1:1) solution to produce a 1 mM primary stock solution. From each of these stock solutions, 1 mL was transferred into a volumetric flask and diluted with methanol to reach a final working solution concentration of 5 μM in methanol of IS\_Cer\_16:0-d31, IS\_FA\_17:0-d33, IS\_LPC\_14:0-d42, IS\_PA\_34:1-d31, IS\_PC\_34:1-d31, IS\_PE\_34:1-d31, IS\_PG\_34:1-d31, IS\_PI\_34:1-d31, IS\_PS\_28:0-d54, IS\_SM\_34:1-d31, IS\_TG\_45:0-d87.

ACYL-CARNITINE-IS: the acyl-carnitine stable isotope-labelled internal standard was prepared by dissolving each powdered stock into methanol to achieve a 5 mM stock solution. Taking 1 mL of the IS\_Car\_4:0-d7 and IS\_Car\_16:0-d3 stock solutions and diluting these into methanol until a final working solution of 5 μM was achieved for IS\_Car\_4:0-d7 and IS\_Car\_16:0-d3.

LIPID-QC: the lipid quality control standards were prepared by diluting each lipid mix to achieve a 50 μg/mL working stock solution in propan-2-ol: acetonitrile: water (2:1:1, respectively).

#### *4.2. Extraction*

Lipids were isolated comparing two methods; firstly, a novel protein-precipitation liquid extraction and secondly the liquid–liquid extraction previously described by Folch and colleagues [12] in an adapted version as we described previously [31]. Tissue quantities ranged from ~2–50 mg and fluid samples from 10–50 μL (e.g., plasma/serum) were tested (data not shown here).

#### 4.2.1. Protein Precipitation Liquid Extraction Protocol (PPLE)

The protein-precipitation liquid extraction protocol was as follows: the tissue samples were weighed (NB. fluid samples were pipetted) and transferred into a 2 mL screw cap Eppendorf plastic tube (Eppendorf, Stevenage, UK) along with a single 5 mm stainless steel ball bearing. Immediately, 400 μL of chloroform: methanol (2:1, respectively) solution was added to each sample, followed by thorough mixing. The samples were then homogenised in the chloroform: methanol (2:1, respectively) using a Bioprep 24-1004 homogenizer (Allsheng, Hangzhou, China) run at speed; 4.5 m/s, time; 30 s for 2 cycles. Then, 400 μL of chloroform, 100 μL of the LIPID-IS (5 μM in methanol) and 100 μL of the CARNITINE-IS (5 μM in methanol) was added to each sample. The samples were homogenised again using a Bioprep 24-1004 homogenizer run at speed; 4.5 m/s, time; 30 s for 2 cycles. To ensure fibrous material was diminished, the samples were sonicated for 30 min in a water bath sonicator (Advantage-Lab, Menen, Belgium). Then, 400 μL of acetone was added to each sample. The samples were thoroughly vortexed and centrifuged for 10 min at ~20,000× *g* to pellet any insoluble material at the bottom of the vial. The single layer supernatant was pipetted into separate 2 mL screw cap amber-glass auto-sampler vials (Agilent Technologies, Cheadle, UK); being careful not to break up the solid pellet at the bottom of the tube. The organic extracts (chloroform, methanol, acetone composition, ~1.4 mL) were dried down to dryness using a Concentrator Plus system (Eppendorf, Stevenage, UK) run for 60 min at 60 ◦C. The samples were reconstituted in 100 μL of 2:1:1 (propan-2-ol, acetonitrile and water, respectively) then thoroughly vortex. The reconstituted sample was transferred into a 250 μL low-volume vial insert inside a 2 mL amber glass auto-sample vial ready for liquid chromatography with mass spectrometry detection (LC–MS) analysis.

#### 4.2.2. Folch Liquid–Liquid Extraction Protocol (Folch LLE)

The Folch liquid–liquid extraction protocol is as follows: the tissue samples were weighed (NB. fluid samples were pipetted) and transferred into a 2 mL screw cap Eppendorf plastic tube (Eppendorf, Stevenage, UK) along with a single 5 mm stainless steel ball bearing. Immediately, 400 μL of chloroform: methanol (2:1, respectively) solution was added to each sample, followed by thorough mixing. The samples were then homogenised in the chloroform: methanol (2:1, respectively) using a Bioprep 24-1004 homogenizer (Allsheng, Hangzhou, China) run at speed; 4.5 m/s, time; 30 s for 2 cycles. Then, 400 μL of chloroform, 100 μL of the LIPID-IS (5 μM in methanol) and 100 μL of the ACYL-CARNITINE-IS (5 μM in methanol) was added to each sample. The samples were homogenised again using a Bioprep 24-1004 homogenizer run at speed; 4.5 m/s, time; 30 s for 2 cycles. To ensure fibrous material was diminished, the samples were sonicated for 30 min in a water bath sonicator. Then, 400 μL of HPLC water was added to each samples. The samples were thoroughly vortexed and centrifuged for 10 min at ~20,000 g to separate the two immiscible fractions. The organic fractions (the lower layer, mostly chloroform; ~700 μL) and aqueous fractions (the upper layer, methanol and water; ~700 μL) were pipetted into separate 2 mL screw cap amber-glass auto-sampler vials (Agilent Technologies, Cheadle, UK); being careful not to break up the solid pellet between the layers. To ensure complete lipid isolation a double extraction protocol was followed; 1 mL of chloroform: methanol (2:1, respectively) solution was added to each sample, along with 400 μL of HPLC water. The samples were thoroughly vortexed and centrifuged for 10 min at ~20,000× *g*. The organic fractions and aqueous fractions were pipetted into the corresponding 2 mL screw cap amber-glass auto-sampler vials containing the initial extracts (again being careful not to break up the solid pellet between the layers). The combined organic extracts (~1.4 mL) were dried down to dryness using a Concentrator Plus system (Eppendorf, Stevenage, UK) run for 60 min at 60 ◦C. The samples were reconstituted in 100 μL of 2:1:1 (propan-2-ol, acetonitrile and water, respectively) then thoroughly vortex. The reconstituted sample

was transferred into a 250 μL low-volume vial insert inside a 2 mL amber glass auto-sample vial ready for liquid chromatography with mass spectrometry detection (LC–MS) lipidomics analysis.

#### *4.3. LC–MS Method*

Full chromatographic separation of intact lipids was achieved using a Shimadzu HPLC System (Shimadzu UK Limited, Milton Keynes, UK) with the injection of 10 μL onto a Waters Acquity UPLC® CSH C18 column (Waters, Hertfordshire, UK); 1.7 μm, I.D. 2.1 mm × 50 mm, maintained at 55 ◦C. Mobile phase A was 6:4, acetonitrile and water with 10 mM ammonium formate. Mobile phase B was 9:1, propan-2-ol and acetonitrile with 10 mM ammonium formate. The flow was maintained at 500 μL per minute through the following gradient: 0.00 min\_40% mobile phase B; 0.40 min\_43% mobile phase B; 0.45 min\_50% mobile phase B; 2.40 min\_54% mobile phase B; 2.45 min\_70% mobile phase B; 7.00 min\_99% mobile phase B; 8.00 min\_99% mobile phase B; 8.3 min\_40% mobile phase B; 10 min\_40% mobile phase B. The sample injection needle was washed using 9:1, 2-propan-2-ol and acetonitrile. The mass spectrometer used was the Thermo Scientific Exactive Orbitrap with a heated electrospray ionisation source (Thermo Fisher Scientific, Hemel Hempstead, UK). The mass spectrometer was calibrated immediately before sample analysis using positive and negative ionisation calibration solution (recommended by Thermo Scientific). Additionally, the heated electrospray ionisation source was optimised at 50:50 mobile phase A to mobile phase B for spray stability (capillary temperature; 300 ◦C, source heater temperature; 420 ◦C, sheath gas flow; 40 (arbitrary), auxiliary gas flow; 15 (arbitrary), spare gas; 3 (arbitrary), source voltage; 4 kV. The mass spectrometer scan rate set at 4 Hz, giving a resolution of 25,000 (at 200 *m*/*z*) with a full-scan range of *m*/*z* 100 to 1800 with continuous switching between positive and negative mode.

#### *4.4. Data Processing*

Thermo Xcalibur Quan Browser (Thermo Fisher Scientific, Hemel Hempstead, UK) data processing involved the integration of the internal standard extracted ion chromatogram (EIC) peaks at the expected retention times (see Table 2). The EIC were selected from the ionisation mode for each analyte class; the ionisation mode is dependent on the molecular chemistry of the analytes, i.e., basic chemical groups ordinarily result in positive ionisation (e.g., [M+H]+, M+H-H2O]<sup>+</sup>, [M+Na]+, [M+NH4] <sup>+</sup>, [M+K]+) whereas acidic chemical groups typically result in negative ionisation (e.g., [M-H]−).

**Table 2.** This table shows the stable isotope-labelled internal standards with their ionisation products (i.e., [M+H]+, M+H-H2O]<sup>+</sup>, [M+Na]+, [M+NH4] <sup>+</sup>, [M+K]+, [M-H]−) and primary ionisation mode (positive; +ve or negative; −ve), along with their retention time (minutes). Butyryl-d7-L-carnitine (abbreviated to IS\_Car\_4:0-d7), N-tetradecylphosphocholine-d42 (abbreviated to IS\_LPC\_14:0-d42), hexadecanoyl-L-carnitine-d3 (abbreviated to IS\_Car\_16:0-d3), heptadecanoic-d33 acid (abbreviated to IS\_FA\_17:0-d33), 1,2-dimyristoyl-d54-sn-glycero-3-[phospho-L-serine] (abbreviated to IS\_PS\_28:0-d54), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphoinositol (abbreviated to IS\_PI\_34:1-d31), N-palmitoyl-d31-D-erythro-sphingosylphosphorylcholine (abbreviated to IS\_SM\_34:1-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-[phospho-rac-(1-glycerol)] (abbreviated to IS\_PG\_34:1-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphate (abbreviated to IS\_PA\_34:1-d31), N-palmitoyl-d31-D-erythro-sphingosine (abbreviated to IS\_Cer\_16:0-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphocholine (abbreviated to IS\_PC\_34:1-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphoethanolamine (abbreviated to IS\_PE\_34:1-d31), glyceryl tri(pentadecanoate-d29) (abbreviated to IS\_TG\_45:0-d87).



**Table 2.** *Cont*.

As shown in the table above (see Table 2), the internal standards have multiple ionisation products, these are the result of numerous ionisation mechanism (for example IS\_TG\_45:0-d87 having different adducts: [M+H]+, [M+Na]+, [M+K]<sup>+</sup> and [M+NH4] <sup>+</sup>, present) as well as an isotopic distribution (e.g., IS\_TG\_45:0-d87 having either the expected eighty-seven or fewer deuterium atoms present) all reasonably expected ions were included into the EIC for each internal standard (see Figure 4).

**Figure 4.** This figure shows a spiked (5 μM in methanol) commercial human plasma extracted ion chromatogram (EIC) for the stable isotope-labelled internal standards (lipids and acyl-carnitines): butyryl-d7-L-carnitine (abbreviated to IS\_Car\_4:0-d7); area ~5.9 <sup>×</sup> 106 counts, N-tetradecylphosphocholine-d42 (abbreviated to IS\_LPC\_14:0-d42); are ~2.1 <sup>×</sup> 108 counts, hexadecanoyl-L-carnitine-d3 (abbreviated to IS\_Car\_16:0-d3); area ~5.8 <sup>×</sup> 108 counts, heptadecanoic-d33 acid (abbreviated to IS\_FA\_17:0-d33); area ~1.1 <sup>×</sup> 104 counts, 1,2-dimyristoyl-d54 sn-glycero-3-[phospho-L-serine] (abbreviated to IS\_PS\_28:0-d54); area ~2.2 <sup>×</sup> 107 counts, 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphoinositol (abbreviated to IS\_PI\_34:1-d31): area ~1.1 <sup>×</sup> 107 counts, N-palmitoyl-d31-D-erythro-sphingosylphosphorylcholine (abbreviated to IS\_SM\_34:1-d31); 1.4 <sup>×</sup> 108 counts, 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-[phospho-rac-(1-glycerol)] (abbreviated to IS\_PG\_34:1-d31); area ~6.1 <sup>×</sup> 107 counts, 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphate (abbreviated to IS\_PA\_34:1-d31); area ~1.3 <sup>×</sup> 107 counts, N-palmitoyl-d31-D-erythro-sphingosine (abbreviated to IS\_Cer\_16:0-d31); area ~2.2 <sup>×</sup> 108 counts, 1-palmitoyl-d31-2-oleoyl-sn-glycero-3 phosphocholine (abbreviated to IS\_PC\_34:1-d31); area ~2.5 <sup>×</sup> 108 counts, 1-palmitoyl-d31-2-oleoylsn-glycero-3-phosphoethanolamine (abbreviated to IS\_PE\_34:1-d31); area ~6.3 <sup>×</sup> 107 counts, glyceryl tri(pentadecanoate-d29) (abbreviated to IS\_TG\_45:0-d87); area ~3.2 <sup>×</sup> 107 counts.

The adduct composition of the total EIC produced from each of the ionisation mechanisms are shown in the figure below (see Figure 5).

**Figure 5.** This figure shows the intensity of the extracted ion chromatogram for each stable isotope-labelled internal standard along with the ionisation adduct composition: [M+H]+, [M+H-H2O]<sup>+</sup>, [M+NH4] <sup>+</sup>, [M+Na]+, [M+K]<sup>+</sup> and [M-H]−. Butyryl-d7-L-carnitine (abbreviated to IS\_Car\_4:0-d7), N-tetradecylphosphocholine-d42 (abbreviated to IS\_LPC\_14:0-d42), hexadecanoyl-L-carnitine-d3 (abbreviated to IS\_Car\_16:0-d3), heptadecanoic-d33 acid (abbreviated to IS\_FA\_17:0-d33), 1,2-dimyristoyl-d54-sn-glycero-3-[phospho-L-serine] (abbreviated to IS\_PS\_28:0-d54), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphoinositol (abbreviated to IS\_PI\_34:1-d31), N-palmitoyl-d31-D-erythro-sphingosylphosphorylcholine (abbreviated to IS\_SM\_34:1-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-[phospho-rac-(1-glycerol)] (abbreviated to IS\_PG\_34:1-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphate (abbreviated to IS\_PA\_34:1-d31), N-palmitoyl-d31-D-erythro-sphingosine (abbreviated to IS\_Cer\_16:0-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphocholine (abbreviated to IS\_PC\_34:1-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphoethanolamine (abbreviated to IS\_PE\_34:1-d31), glyceryl tri(pentadecanoate-d29) (abbreviated to IS\_TG\_45:0-d87).

The data processing also involved the integration of the individual lipid (and derivatives) species at their expected retention time (see Supplementary Table S2) allowing for a maximum of ±0.1 min of retention time drift: any retention time drift greater than ±0.1 min resulted in the exclusion of the analyte leading to a 'Not Found' result (i.e., zero concentration). A list of the analyte classes along with the number of species detected within each class are shown in the table below (see Table 3). The expected adducts for each analyte class and the internal standard used for semi-quantitation are also shown.

The lipid quality control (QC) standards were analysed with each batch of samples, these QC standards were used to check the retention times for the analytes ensuring that isobaric analytes were separated and expected analyte retention times remained robust.

Through the Thermo Xcalibur Quan Browser software, the responses of the analytes were normalised to the relevant internal standard response (producing area ratios) (see Supplementary Table S2), these area ratios corrected the intensity for any extraction and instrument variations. The area ratios were then blank corrected where intensities less than three times the blank samples were set to

a 'Not Found' result (i.e., zero concentration). The accepted area ratios were then multiplied by the concentration of the internal standard to give the analyte concentrations. The results for fluid samples were expressed in molar concentrations (typically μM or nM). For tissue samples, the calculated concentrations of the analytes were then divided by the amount of tissue (in mg) used in the extraction protocol to give the final results in μM per mg of tissue extracted (μM/mg).

**Table 3.** This table shows the lipid classes detected with this LC–MS lipidomics method. The number of species per lipid class and the measured adducts (protonated: [M+H]+, deprotonated: [M-H]−, protonated with water loss: [M+H-H2O]<sup>+</sup>, sodiated: [M+Na]+, potasiated: [M+K]+, ammoniated: [M+NH4] <sup>+</sup>) are also shown. The internal standard used for semi-quantification are also shown: butyryl-d7-L-carnitine (abbreviated to IS\_Car\_4:0-d7), N-tetradecylphosphocholine-d42 (abbreviated to IS\_LPC\_14:0-d42), hexadecanoyl-L-carnitine-d3 (abbreviated to IS\_Car\_16:0-d3), heptadecanoic-d33 acid (abbreviated to IS\_FA\_17:0-d33), 1,2-dimyristoyl-d54-sn-glycero-3-[phospho-L-serine] (abbreviated to IS\_PS\_28:0-d54), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphoinositol (abbreviated to IS\_PI\_34:1-d31), N-palmitoyl-d31-D-erythro-sphingosylphosphorylcholine (abbreviated to IS\_SM\_34:1-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-[phospho-rac-(1-glycerol)] (abbreviated to IS\_PG\_34:1-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphate (abbreviated to IS\_PA\_34:1-d31), N-palmitoyl-d31-D-erythro-sphingosine (abbreviated to IS\_Cer\_16:0-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphocholine (abbreviated to IS\_PC\_34:1-d31), 1-palmitoyl-d31-2-oleoyl-sn-glycero-3-phosphoethanolamine (abbreviated to IS\_PE\_34:1-d31), glyceryl tri(pentadecanoate-d29) (abbreviated to IS\_TG\_45:0-d87).


#### *4.5. Animal Intervention*

Sprague–Dawley rats (Harlan, IN, USA) were overfed using one of four experimental diets (n = 6–9 per group) at 17% above matched growth via an intragastric cannula surgically inserted as previously described [26]. Animals had ad libitum access to water throughout the experiments. The four experimental diets were 70% fat (% energy) including different amounts of medium chain triacylglycerides oil (MCT), beef tallow and corn oil; the fat composition of each diet are shown in the table below (see Table 4).

Protein (19% whey protein), vitamin and mineral contents were the same in all diets. Diets were formulated to meet the caloric and nutritional recommendations established by the National Research Council (NRC), but were fed at a level that exceeded the recommended caloric intake by 17% to increase weight gain and adiposity and produce steatohepatitis.

Liver tissue was collected after 21 days. All experimental procedures were ethically approved by the Institutional Animal Care and Use Committee at the University of Arkansas for Medical Science.


**Table 4.** This table shows the dietary fat composition of each of the four experimental diets fed to Sprague–Dawley rats (n = 8–9 per group). MCT: medium chain triglyceride oil.

#### **5. Conclusions**

This lipidomics protocol has been developed to quantify lipids across a broad range of hydrophobicities, from acyl-carnitines through to long chain glycerolipids. The extraction method produces a single liquid supernatant phase ideal for high-throughput workflows with an increased extraction capability over the frequently published liquid–liquid extraction previously published by Folch and colleagues [12].

Following the establishment and validation of this method, we applied it to a ruminant fat dose response dietary intervention in Sprague–Dawley rats, where we found 100 lipid species correlated strongly with the composition of ruminant fat within the diet.

It has been previously suggested that dietary ruminant fat is beneficial/protective in type 2 diabetes [32], the results presented in this manuscript suggest possible target mechanisms that need to be examined could include ceramide fatty acid compositions, cardiolipin remodeling, sphingomyelins and/or triacylglycerides concentration (particularly unsaturated odd chain species) and their associated fatty acid compositions, as well as the liver total lipid content.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2218-1989/10/7/296/s1, Figure S1: This figure shows the comparison between the two lipid extraction techniques regarding their extraction efficiency on each lipid class detected in the rat liver samples (Folch liquid–liquid extraction with a compositions of chloroform: methanol: water, ~7:3:4, and Protein precipitation liquid extraction with a composition of chloroform: methanol: acetone, ~7:3:4). n = 34 rat liver samples per extraction method. The intensity of the lipids were measured by liquid chromatography with mass spectrometry. The significance of the difference between the two extraction protocols are shown by the *p*-value star system; where *p* ≤ 0.05 was considered statistically significant (\* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001). Error bars represent ± standard deviation. Table S1: This table shows the comparison between the two lipid extraction techniques regarding their extraction efficiency on the total intensity of each lipid class detected in the rat liver samples by each sample extraction method (Folch-LLE: Folch liquid–liquid extraction with a compositions of chloroform: methanol: water, ~7:3:4, and PPLE: protein precipitation liquid extraction with a composition of chloroform: methanol: acetone, ~7:3:4). n = 34 rat liver samples per extraction method. The percentage increase the PPLE method is over the Folch-LLE method is shown (% diff.) along with the *p*-value resulting from a *t*-test (*p* < 0.05 designates statistical significance are in bold & shaded). The total number of lipid species detected and pass the quality control process are also shown. Table S2: This table shows the lipids quantified in this LC–MS method, along with the ionisation mode (either positive; +ve, or negative; −ve), the detected ion (*m*/*z*), the expected retention time (minutes) and the internal standard used for normalisation and quantification. Lipid are shown in their shorthand notations with the number of carbons and unsaturated bonds in the fatty acid moiety separated by a colon; acyl-carnitines (Carn), ceramides (Cer), cardiolipins (CL), diacylglycerols (DG), gangliosides (GM1), hexosylceramides (Hex-Cer), lyso-phosphatidylcholines (LPC), lyso-phosphatidyethanolamines (LPE), lyso-phosphatidylinositols (LPI), lyso-cardiolipins (Lyso\_CL), phosphatidic acids (PA), phosphatidylcholines (PC), phosphatidylethanolamines (PE), phosphatidylglycerol (PG), phosphatidylinositols (PI), phosphatidylserines (PS), sulfatides (S), sphingomyelins (SM), triacylglycerides (TG).

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

**Funding:** This research was funded by BBSRC, grant number BB/M027252/1, USDA, grant number ACNC-USDA-CRIS 6251-51000-005-03S, and COST Action, grant number CM0603.

**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*
