Nonalcoholic fatty liver disease stratification by liver lipidomics.

Nonalcoholic fatty liver disease (NAFLD) is a common metabolic dysfunction leading to hepatic steatosis. However, NAFLD's global impact on the liver lipidome is poorly understood. Using high-resolution shotgun mass spectrometry, we quantified the molar abundance of 316 species from 22 major lipid classes in liver biopsies of 365 patients, including nonsteatotic patients with normal or excessive weight, patients diagnosed with NAFL (nonalcoholic fatty liver) or NASH (nonalcoholic steatohepatitis), and patients bearing common mutations of NAFLD-related protein factors. We confirmed the progressive accumulation of di- and triacylglycerols and cholesteryl esters in the liver of NAFL and NASH patients, while the bulk composition of glycerophospho- and sphingolipids remained unchanged. Further stratification by biclustering analysis identified sphingomyelin species comprising n24:2 fatty acid moieties as membrane lipid markers of NAFLD. Normalized relative abundance of sphingomyelins SM 43:3;2 and SM 43:1;2 containing n24:2 and n24:0 fatty acid moieties, respectively, showed opposite trends during NAFLD progression and distinguished NAFL and NASH lipidomes from the lipidome of nonsteatotic livers. Together with several glycerophospholipids containing a C22:6 fatty acid moiety, these lipids serve as markers of early and advanced stages of NAFL.

Nonalcoholic fatty liver disease (NAFLD) is a metabolic dysfunction histologically characterized by hepatic fat accumulation (hepatic steatosis) in the absence of heavy alcohol consumption in past medical history (1). NAFLD affects up to 30% of adults and up to 80% of obese and diabetic individuals worldwide (2). The prevalence and severity of NAFLD are higher in men although in postmenopausal women, the NAFLD rate increases (3). NAFLD is subdivided into nonalcoholic fatty liver (NAFL), nonalcoholic steatohepatitis (NASH), cirrhosis, and hepatocellular carcinoma (4, 5). In contrast to NASH, hepatic steatosis in NAFL may occur with no significant inflammation.

The molecular background and pathophysiology of NAFL and why and how it progresses to NASH are poorly understood (6, 7, 8, 9). While the intracellular accumulation of triacylglycerols (TG), diacylglycerols (DG) and free fatty acids (FFA) is a metabolic hallmark of NAFLD, it is unclear whether the disease also alters a broader scope of lipids (10, 11, 12, 13, 14, 15, 16). Mutations in NAFLD risk factors, e.g., PNPLA3 or MBOAT7, affect the liver lipidome (17, 18, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31). It is also conceivable that lipidome remodeling may contribute to or be associated with the onset and propagation of NAFLD (15). However, lipidomic evidence is mostly semiquantitative and based on limited lipid class coverage (32). Previous studies were often focused on the composition of energy storage lipids (12, 13, 33), while membrane and signaling lipids received less attention (28, 34, 35). Chiappini et al. (36) used TOF-SIMS imaging to analyze 104 lipid species in biopsies of 61 NAFLD patients and suggested a NASH lipidomic signature comprising 32 lipid species, although the analysis was biased to most abundant and best ionized classes, e.g., phosphatidylcholines (PC) and TG. The study by Gorden et al. (13) covered 186 lipids in liver biopsies and plasma of 91 patients and revealed a combination of plasma biomarkers (including polyunsaturated glycero- and glycerophospholipids together with long-chain ceramides (Cer)) that distinguished NASH from NAFL. However, lipid markers and their fold changes between disease states reported by independent studies were not concordant. Also, if NAFLD globally alters the liver lipidome and if changes in the abundance of glycerophospho- and sphingolipids (for convenience, here we termed them as membrane lipids) corroborate major clinical indices and correlate with disease severity is an open question, particularly because no reference values for individual liver lipids and their physiological variation are available.

To better understand how NAFLD is transforming the human liver lipidome, we assembled a cohort of 365 histologically characterized biopsies reflecting its progression from nonsteatotic obesity to overt NASH, together with appropriate nonsteatotic and nonobese controls. We then used shotgun mass spectrometry (37, 38) to systematically quantify the molar abundance of 316 species from 22 major lipid classes that encompassed membrane and energy storage lipids including cholesterol. Biclustering analysis of the curated lipidomics dataset recognized several signatures comprising specific membrane lipids that enabled patient's stratification at different stages of NAFLD independently of progressive accumulation of DG and TG.

Materials and Methods

Cohort recruitment, study design, and ethic approval

The study protocol accords the ethical guidelines of the 1975 Declaration of Helsinki and was approved by the authority of Universität Kiel (D425/07, A111/99) before the study commenced. All patients had given their written informed consent. In total, 365 individuals (124 males; 241 females) from 17 to 85 years of age and whose BMI was in the range of 14.8–83.6 (kg/m2) were recruited within the time period of 2007–2016. NASH and NAFL were defined by the NAFLD activity score (NAS) as described (39). Phenotyping of the entire cohort was performed using standardized histology protocol (17) in a blinded fashion by a board-certified surgical pathologist (C. R.) having the specialization in hepatopathology (details are in supplemental data, Histology).

Alcohol consumption was assessed by self-reporting; subjects with average alcohol consumption of more than 30 g/day in men or 20 g/day in women (an equivalent of three and two standard alcoholic drinks per day, respectively) were not enrolled (40). The collected metadata included age, sex, BMI, blood test, including gamma-glutamyl transferase (GGT), medication taken by each patient, and mutation status of the following genes: PNPLA3, TM6SF2, MBOAT7, HSD17B13, SERPINA1 (SERPINA S and SERPINA Z). Evidence for insulin resistance was not available.

Common chemicals and lipid standards

Synthetic lipid standards (see supplemental data, Common chemicals and lipid standards for complete list) were purchased from Avanti Polar Lipids (Alabaster, AL, USA). Individual standards were mixed and diluted with methyl-tert-butyl ether (MTBE)/methanol (MeOH) 10:3 (v/v) (see supplemental data, Common chemicals and lipid standards for details).

Annotation of lipid classes and species

Glycerolipids are referred to TG and DG; glycerophospholipids and lyso-glycerophospholipids to phosphatidic acids (PA), phosphatidylinositols (PI), phosphatidylserines (PS), phosphatidylglycerols (PG), phosphatidylethanolamines (PE), phosphatidylcholines (PC), ether phosphatidylethanolamines (PE O-), ether phosphatidylcholines (PC O-), lyso-phosphatidic acids (LPA), lyso-phosphatidylinositols (LPI), lyso-phosphatidylcholines (LPC), and lysophosphatidylethanolamines (LPE); sphingolipids to ceramides (Cer) and sphingomyelins (SM); sterols to cholesterol (Chol) and cholesteryl esters (CE). Species of glycero- and glycerophospholipids and cholesteryl esters are annotated as : in both (or, for lyso-lipids and cholesteryl esters, in one) fatty acid or fatty alcohol moieties (moiety). Sphingolipids are annotated as : ; at the ceramide backbone.

Sample preparation for shotgun lipidomics

Biopsies (wet weight of 4.2–21.9 mg) were shock-frozen in liquid nitrogen ensuring an ex vivo time of less than 40 s and stored at –80°C freezer. Prior to lipid extraction, the tissues were homogenized in 300 μl of isopropanol using zirconium beads; the total protein content was determined by Pierce 660 assay (Thermo Fisher Scientific, USA). Lipids were extracted from aliquots containing an equivalent of 50 μg of total protein by adding 700 μl of MTBE/MeOH 10:3 (v/v) containing the internal standard mix (41, 42) (see details in supplemental data, Sample preparation). After evaporation of the organic phase, lipid extracts were reconstituted in 600 μl of 2:1 (v/v) MeOH / CHCl3 and stored at –20°C. Ten microliters of a lipid extract were diluted with 90 μl of the spray solution (4:2:1 isopropanol/MeOH/CHCl3 (v/v/v) containing 7.5 mM ammonium formate) for mass spectrometric analysis. Samples were analyzed in technical duplicates.

Lipid identification and quantification by shotgun mass spectrometry

The mass spectrometric analysis was performed on a Q Exactive instrument (Thermo Fisher Scientific, Bremen, Germany) equipped with a robotic nanoflow ion source TriVersa NanoMate (Advion BioSciences, Ithaca, NY) using nanoelectrospray chips with spraying nozzle diameter of 4.1 μm. The ion source was controlled by the Chipsoft 8.3.1 software (Advion BioSciences). Ionization voltage was +0.96 kV in positive and −0.96 kV in negative mode; backpressure was 1.25 psi in both modes (43). Temperature of the ion transfer capillary was 200 °C; S-lens RF level was 50%. FT MS spectra were acquired within the range of m/z 400–1,000 in positive and m/z 350–1,000 in negative ion mode at the target mass resolution of R m/z 200=140,000; automated gain control (AGC) of 3 × 106 and maximal injection time of 3 s. In both modes t-SIM spectra were acquired at the same mass resolution and m/z range as above; AGC of 5 × 104; maximum injection time of 650 ms; width of isolation window of 20 Th. The inclusion list of masses targeted by t-SIM started at m/z 355 in negative and m/z 405 in positive ion mode and other masses were computed by adding 10 Da increment (i.e., m/z 355, 365, 375) up to m/z 1,005. Free cholesterol was quantified by parallel reaction monitoring (PRM) FT MS/MS (41, 44, 45) during the same analysis. The number of micro-scans was set to 1; width of precursor isolation window of 0.8 Th; normalized collision energy (nCE): 12.5%; AGC: 5 × 104 and maximum injection time of 3 s.

Raw t-SIM spectra were subjected to repetition rate filtering by PeakStrainer software (46) and then stitched together by SIMStitcher software (47). Lipids were identified by LipidXplorer software (48) by accurately determined m/z (mass accuracy better than 5 ppm) and quantified by comparing the isotopically corrected abundances of their molecular ions with the abundances of internal standards of the same lipid class. All internal standards were detected in all samples. MS2 validation of selected species by HCD FT MS/MS was applied as described in supplemental data, MS2 validation.

Raw data processing

The LipidXplorer output was processed using several steps. Technical replicates were averaged and spectra of QC samples (details in supplemental data, Pilot sub-cohort and Quality Control) set aside. In each sample lipid abundances for which the standard deviation (SD) exceeded 40% were set to Not-a-Number (NaN). Lipids whose abundance exceed the minimal value determined for this species by less than 2-fold were exempted from SD filtering (supplemental data, Raw data processing and supplemental Fig. S1). Also, cholesterol values were not SD-filtered because they were determined by PRM. Next, we grouped the patient samples according to their disease status (normal control, healthy obese, NAFL, NASH, and none). The abundance of lipid species that were detected in less than 15% of all patient samples of the same group was set to zero. Finally, we applied plate bias correction using the ComBat approach (49) as detailed in supplemental data, Raw data processing. The final dataset comprised the molar abundance of 316 lipid species from 22 lipid classes in 365 patients (supplemental dataset S1).

Bioinformatics and data mining

Biclustering analysis was performed using the isa (50) and qubic (51) algorithms. To eliminate redundancy of biclusters, we screened them using a consensus approach (see supplemental data, Bioinformatics and data mining).

Classifications were performed using random forest classifiers from the “randomForest” R package. Hundred classifiers were trained per scenario on bootstrapped lipidomics data. The models were evaluated using precision recall (PR) and receiver operator characteristics (ROC) area under the curve (AUC) values on the test set. Feature importance was computed using gini index averaged over all bootstrapped models (For further details, see supplemental data, Bioinformatics and data mining).

Lipidomics impact of mutations

The impact of mutations was probed within each group (supplemental data, Genotyping), in which the lipidomes of homozygotes, heterozygotes, and normal genotypes were compared pairwise (details are in supplemental data, Statistical analysis of mutation-correlated lipidome changes).

Results

Study cohort

Percutaneous or surgical liver biopsies were taken from in total of 365 Caucasian patients (124 males and 241 females) admitted to the University Hospital Schleswig-Holstein, Kiel, Germany during 2007–2016. Control biopsies were collected during gastrointestinal surgery for pathologies having no direct association with NAFLD. Each biopsy was examined by the same surgical pathologist (C. R.) (17) and assessed according to the NAFLD activity score (NAS) (39).

Patients were further sorted into four basal and one additional group (52) as outlined below. Biopsies were classified according to the presence of steatosis and inflammation as the main criteria to differentiate simple steatosis (NAFL) from steatohepatitis (NASH) (52). Therefore, individuals whose biopsies were having histological fat content ≤5%; BMI <30 kg/m2; no histologically proven liver pathology; no lobular inflammation and absence of significant (stage 2 or higher) fibrosis formed the normal control (NC) group. Healthy obese (HO) group comprised patients having no NAFLD per clinical and histological evidence, yet their BMI exceeded 30 kg/m2. Patients were assigned to NASH and NAFL groups if the histological examination revealed steatosis with or without inflammation, respectively. In total, 337 patients were sorted into NC (n = 49), HO (n = 51), NAFL (n = 143), and NASH (n = 94) groups (Table 1; supplemental Table S1, age and BMI distribution is in supplemental Fig. S4).

Table 1

Clinicopathological characteristics of the study cohort

Group	Number of patients			NAS Parametersa				Fibrosisa	Mean age, years	Mean BMI kg/m2
	Total	Female	Male	NAS Fat	NAS Ballooning	NAS Inflammation	NAS Average
NC	49	25	24	0	0	0	0	0	67.7	24.1
HO	51	47	4	0	0	0	0	0	42.3	46.7
NAFL	143	94	49	1	0	0	0	0.35	45.6	47.1
NASH	94	62	32	2	0	1	4	1	45.1	49.3
None	28	13	15	n.a.	n.a.	n.a.	n.a.	n.a.	57.1	33.7

Median values, in arbitrary units assigned according to (39).

Within the cohort, 28 patients who did not receive NAS because of serious liver conditions, such as histologically confirmed cancer with abnormal and higher than control levels of GGT, alkaline phosphatase (AP), and bilirubin, were assigned to the additional group “None.” Their lipidomes helped to delineate the impact of a broader spectrum of liver diseases different from NAFLD. The impact of obesity (53) was evaluated by comparing lipidomes of NC and HO groups. Since mean BMI of HO, NAFL, and NASH groups was similar, their lipidomes could be compared with no further adjustment. The mean age of NC group was 12 years higher than of NAFL/NASH/None groups because of fewer younger people undergoing gastrointestinal surgery.

Shotgun lipidomics of liver biopsies

Although surgical biopsies were similar in size, they varied in weight by as much as 5-fold and also differed in fat and fibrotic content (Table 1). Since they were unique and histologically inhomogeneous, each specimen was analyzed as a single sample with no independently processed replicates. In an aliquot of each biopsy lysate, we first determined the total protein content and used it for subsequent normalization of molar lipid abundances.

To quantify the lipidomes, we employed high-resolution shotgun Fourier transform mass spectrometry (FT MS) that relies on direct nanoflow infusion of total lipid extracts (37, 38). During shotgun analysis, internal standards spiked into biopsy lysates prior lipid extraction are ionized together with endogenous lipids and enable their absolute (molar) quantification (17, 42, 54, 55, 56). Because of high and variable content of TG together with abundant chemical background, we analyzed liver extracts by the method of t-SIM (for targeted single ion monitoring) (47). In each analysis FT MS spectra were successively acquired from partially overlapping m/z windows of 20 Th and then these spectra were stitched together into a master spectrum by SIMStitcher software. Although typical t-SIM acquisition (supplemental Fig. S5) required ca. 11 min per sample (compared with less than 1 min required for conventional FT MS analysis), it increased the number of quantifiable lipid species by ca. 45% (57).

To ensure spectra acquisition consistency and enable batch correction of lipid abundances, we created a QC standard by pooling equal aliquots of 36 extracts that, according to a preliminary experiment, reflected extreme values of the total lipid content and represented both genders, stages of inflammation and fibrosis. In each aliquot of the QC sample, we quantified 255 lipid species from 19 lipid classes. The molar abundance of more than 98% of lipid species differed by less than 5% and only cholesterol (quantified by the method of PRM (57, 58, 59)) and a few TG species were detected with higher SD (supplemental Fig. S2).

The final dataset covering 365 biopsies comprised normalized molar abundances (in pmol per μg of total protein) of 316 lipid species from 22 lipid classes: Out of the total of 316 lipid species, only 114 were detected in more than 98% of biopsies (supplemental dataset S2).

Liver lipidome of the normal control group

The NC group (n = 49) combined patients with no apparent liver pathology and obesity, whose lipidome we assumed as basal (Fig. 1). The mol% of lipid classes generally corroborates previous reports (13) (supplemental Table S2), although we detected considerably more Cer, PS, SM, and lyso-lipids.

Figure 1

Lipid class composition (in pmol per μg of total protein; log2 scaled ) of liver biopsies in the four main groups of patients. The total abundance of each lipid class was calculated by adding up molar abundances of lipid species. Boxing highlights values between 25% and 75% quartiles; vertical lines connect minimum and maximum values excluding outliers. Filled circles stand for the lipid class abundance in individual biopsies (average of three technical replicates); black lines within boxes indicate median values. Color coding is shown in the inset at the right-hand side.

Since the liver is a hub of lipoprotein biosynthesis, we further examined if the NC lipidome differs from plasma lipidome (42, 56, 60). While both lipidomes have similar mol% proportions of PC, SM, and Chol, plasma contained 10-fold more CE (supplemental Fig. S6). Liver lipidome contains a larger variety of lyso-species (liver-specific are LPG, LPI, LPS, LPA), while LPC is more abundant in mol%. TG and DG are also more abundant in the liver, as well as PE, PG, PI, PS, and Cer. Lipid class profiles in the liver showed no pronounced gender bias (supplemental Fig. S6) that was apparent in healthy plasma (42).

Mutation-associated changes of the liver lipidome

We further evaluated how common mutations in NAFLD risk genes: TM6SF2 variant rs58542926 (hcv89463510), PNPLA3 variant rs738409 (hcv7241), MBOAT7 variant rs641738 (hcv8716820), SERPINA1 variant PiZ (Glu342Lys) rs28929474 (hcv34508510), SERPINA1 variant PiS affected the liver lipidome. Within NC and NASH groups liver lipidomes of homozygotes, heterozygotes and normal genotype carriers were compared pairwise (Table 2). Note that, for consistency, we compared absolute (in moles per μg of total protein) rather than relative (in mol% within each class) lipid abundances, although the latter could reveal more affected species (17) (supplemental Dataset S1).

Table 2

PNPLA3 and MBOAT7 mutation-dependent changes of lipidome in NASH and NC groups

Lipida	Mutation	Group	Mutation Status	Abundance Ratio
CE 22:2	PNPLA3 rs738409	NASH	Heterozygote/No mutation	4.0b
TG 56:7	PNPLA3 rs738409	NASH	Heterozygote/No mutation	6.8c
TG 56:7	PNPLA3 rs738409	NASH	Homozygote/No mutation	15.1c
TG 56:7	PNPLA3 rs738409	NASH	Homozygote/Heterozygote	2.2b
PI 36:1	MBOAT7 rs641738	NC	Heterozygote/No mutation	2.2b
PI 36:1	MBOAT7 rs641738	NC	Homozygote/No mutation	4.5d

Significantly changed lipid in the two compared groups.

Significance P < 0.05.

Significance P < 0.005.

Significance P < 0.01.

MBOAT7 mutation-dependent changes were only detected in the NC group and confined to PI 36:1 (Table 2), consistent with its phosphatidylinositol O-acyltransferase activity and the role in NAFLD pathogenesis (17). In the NASH group PNPLA3 rs738409 mutation elevated the levels of CE 22:2 and TG 56:7, while the total abundance of CE and TG classes did not change (supplemental Figs. S7, S8A,B) (19, 20, 29, 31, 61, 62, 63).

Condition-associated changes in lipid classes

There was no marked difference between the total abundance of lipid classes in disease groups compared with NC (Fig. 1), except significantly higher levels of TG and, to lesser extent, of DG and CE in NAFL and NASH patients (supplemental Figs. S9–S11). Although the mean BMI of HO patients was almost 2-fold higher compared with NC and was as high as in NAFL and NASH patients, the abundance of TG in HO and NC was similar.

In principal component analysis (PCA) plots of patients' lipidomes, the gradient across PC1 (principal component 1) from NC (at the left) to NASH (at the right) (Fig. 2A) reflected the increased abundance of glycerolipids and CE (Fig. 2B). Interestingly, the profiles of hydrocarbon chain length and unsaturation of fatty acid moieties in TG and DG species were the same in all patient groups (Fig. 2C).

Figure 2

Segregation of patient groups by the composition of lipidomes. A: PCA plot for full lipidomes of liver biopsies of 365 patients, whose disease group is indicated by color (coding scheme is in the inset); (B) PCA plot of lipidomes from which glycerolipids and CE were omitted. C: Length and unsaturation of fatty acid moieties in TG and DG species. Circle size and color reflect lipid abundances. D: Similarity network of biclusters. Node size is proportional to the number of patients in the bicluster. Highlighted are four network components (annotated) comprising the largest number of connected biclusters attributed to steatotic (in red and amber for NASH and NAFL) or nonsteatotic patients (in blue and green for NC and HO).

We further looked if TG and DG accumulated in the liver of NAFL and NASH patients were compositionally different from the adipose tissue. The ten most abundant TG species in liver (current dataset) and white adipose tissue (WAT) (64, 65) (supplemental Table S3) were almost the same suggesting that TG metabolism in tissues is not organ-specific. However, they differed from TG in plasma of obese patients (66), presumably because in tissues TG accumulation sequesters lipotoxic FFA, while plasma TG are packed into lipoproteins for transporting via bloodstream (65).

Molecular stratification of patient groups by biclustering

Lipid class composition of liver biopsies (Figs. 1, 2) appears to be conserved and, apart from major increase in glycerolipids and CE, offers limited molecular insight or diagnostic perspective. We reasoned that better molecular stratification of patients could rely on clusters of individual species spanning multiple lipid classes whose molar abundance coherently changes with the disease progression. Specifically, we employed a biclustering ensemble approach (Bioinformatics and data mining), in which lipidome compositions are clustered and the results consolidated into networks of connected components comprising lipid compositions specific for a selection of patients (Fig. 2D; supplemental Figs. S13, S14). Lipid compositions and meta-data could be compared across biclusters and reveal lipid signatures specific for a disease condition or patient group. Within all patients' lipidomes, we recognized two components comprising higher proportion of steatotic (NASH and NAFL) and two components with mostly nonsteatotic (NC and HO) patients (Fig. 2D, supplemental Table S4). Their lipidomes were clearly separated by PCA (supplemental Fig. S12A).

We first tested if lipids identified by biclustering could detail the trend toward steatosis and identify markers distinguishing nonsteatotic from steatotic groups by training random forest-based classifiers. These groups were not separated by common clinical indices, e.g., total bilirubin, GGT, AP, alanine aminotransferase (ALT), and aspartate transaminase (AST) (supplemental Fig. S15B), but were readily distinguished by lipidome compositions (Fig. 3A). Classification performance was evaluated using the precision recall and ROCs (supplemental Fig. S16).

Figure 3

Classification of disease groups by lipid markers. A: Classification of steatotic and nonsteatotic patients reported as an area under the curve (AUC) for precision recall (PR) and receiver operator characteristics (ROC). Classification based on all lipids (upper panel) and only by the ratio of SM 43:1;2 and SM 43:3;2 normalized to SM 38:2;2 (lower panel). B: Changes in the abundance of SM species in steatotic (red background) and nonsteatotic (blue background) components. SM 38:2;2 belongs to none of the two components (gray background). Lines indicate disease progression: NC (blue), HO (green), NAFL (yellow), and NASH (red) patient groups. Axes indicate relative abundances of SM species; data points per species sum up to 100%. C: Medians of the ratios of the abundance of SM 43:1;2 to SM 38:2;2 and SM41:3;2 to SM 38:2;2 for the patient groups, including subgroups of NAFL (D) Classifications of all patient groups against each other using the two SM ratios.

Lipids most significant for the classification were mono- and diunsaturated TG and DG (e.g., TG 50:1, TG 52:2, or DG 34:1) (supplemental Fig. S15A). To further test if unsaturation of fatty acid moieties in TG plays a role, we trained two classifiers with TG having in total 1–4 and 5–8 double bonds, respectively. Both models performed significantly better than random. However, the median Precision Recall Area Under the Curve (PR AUC) of the model with more saturated TG showed better performance (supplemental Fig. S15C, D).

Taken together, our analysis revealed several mono- and diunsaturated DG and TG species as lipid markers of steatosis while common membrane lipids (including glycerophospholipids (GPL), Chol, or Cer) were not in the markers list with one notable exception of two SM species (supplemental Fig. S15A) that also showed clear NAFLD-dependent profile (Fig. 3B).

Sphingomyelin species sharing n24:2 fatty acid moiety are same-class membrane lipid markers of NAFLD

We plotted the relative abundances of all SM species found in steatotic and nonsteatotic groups (Fig. 3B). The abundance of species from the nonsteatotic group (such as SM 42:2;2, SM 40:1;2, SM 42:1;2) either changed marginally or, as for SM 43:1;2, markedly decreased with NAFLD progression with the exception of SM 34:0;2 (Fig. 3B & supplemental Fig. S17).

In contrast, the abundance of SM 41:3;2, SM 43:3;2, and SM 42:4;2 from the steatotic group followed the opposite trend: it consistently increased with NAFLD progression. To elaborate on this finding, we first validated their identification by high-resolution HCD FT MS/MS (supplemental Materials and Methods, MS2 validation). We note that these lipids were previously identified in human plasma by the method of LC-MSn. Major species of SM 41:3;2 and SM 43:3;2 were recognized as SM d17:1/n24:2 and SM d19:1/n24:2, respectively (67); and SM 42:4;2 was SM d18:2/n24:2 (60), suggesting that these three SM may be sharing n24:2 fatty acid moiety. In contrast to SM 41:3;2, another compositionally related yet more saturated sphingomyelin SM 43:1;2 was attributed to nonsteatotic group and its abundance dropped with NAFLD progression.

To validate the molecular composition of liver sphingomyelins, we subjected a few combined extracts of steatotic and nonsteatotic biopsies to targeted LC-MSn (67) (method details are in supplemental data, Identification of SM molecular species by LC-MSn). The analysis (supplemental Table S5) confirmed that SM 34:0;2 belongs to dihydrosphingomyelins, the lipid class associated with fat deposition in organs, including the liver and pancreas (68). SM 42:4;2 and SM 43:1;2 comprised unsaturated (n24:2) and saturated (n24:0) fatty acid moieties, respectively (supplemental Fig. S19; supplemental Table S5). Because of their low abundance, LC-MSn analysis of SM 43:3;2 and SM 41:3;2 was inconclusive.

We noticed that the abundance of SM 41:3;2 and SM 43:3;2 (both comprising n24:2) as compared with SM 41:1;2 (comprising n24:0) followed opposite trends (Fig. 3B), and we examined if their ratio to some unchanged SM species could distinguish different stages of NAFLD. We computed the ratios of SM 43:1;2 (decreasing in NAFLD) and of SM 43:3;2 (increasing in NAFLD) to SM 38:2;2 (SM d18:2/n20:0; see supplemental Table S5) that did not associate with steatotic and nonsteatotic groups and whose abundance was not affected by NAFLD (Fig. 3B). These ratios could be determined directly from a shotgun spectrum of a total lipid extract without prior adjustment to the abundance of internal standards or biopsy size (Fig. 3C). Strikingly, they distinguished nonsteatotic and steatotic subcohorts with only slightly lower specificity than all lipids (Fig. 3A). SM-based classification was unaffected by obesity: it marginally distinguished NC and HO, but readily delineated NC from NAFL and NASH (Fig. 3D). To corroborate these findings, we further examined mol% profiles of species of SM and Cer classes in the four patient groups (supplemental Figs. S17, S18). They confirmed the notable decrease in the marker SM 41:1:2 as well as of the most abundant SM 34:1;2 together with the concomitant increase of SM 41:3;2 and SM 43:3;2 in NAFL and NASH subgroups, as compared with NC and HO. This, however, did not affect the abundance of matching Cer species: the mean value of Cer 34:1;2 was unchanged and none of Cer with odd number of carbon atoms was detectable.

We noted that SM ratios were less specific when we compared HO and NAFL (Fig. 3D). At the same time, HO and NASH were distinguished notably better albeit there was no global difference between their lipidomes except higher abundance of neutral lipids. Therefore, we hypothesized that the cohort of NAFL patients may be compositionally heterogeneous and consist of smaller, yet compositionally distinct subcohorts reflecting some intermediate stages of NAFLD pathogenesis. In this way, SM ratios and likely other lipid markers could reflect the transition from initial and, likely, reversible stage(s) of NAFL towards NASH.

Heterogeneity of NAFL lipidome and its molecular markers

We noticed that a component “Non-steatotic 1” (Fig. 2C) covers the largest number of nonsteatotic (NC+HO) patients, but also some patients with NAFL and NASH. Similarly, “Steatotic 2” mainly covers steatotic (NAFL and NASH) patients and almost no patients of NC and HO groups. We hypothesized that lipidomes of some NAFL patients may have higher similarity to the lipidomes of HO or of NASH. Based on the similarity of lipidome compositions, we used biclustering to divide NAFL patients into four subgroups that clustered together with nonsteatotic (h- for healthy) individuals, with NASH (s- for sick) individuals or neither h- nor s- (hs-) (Fig. 4). To test if this lipidome-based clustering reflects the disease progression also within the NAFL group, we compared the median values of individual histopathological indices (e.g., liver fat mass, ballooning, fibrosis, and inflammation) that were used for calculating NAS. We observed clear disease-related trends for each index and also for the total NAS (supplemental Table S6). Within NAFL subgroups, fat and ballooning expectantly increased from h-NAFL to s-NAFL, whereas fibrosis and inflammation did not, indicating that developing steatosis has not yet led to NASH. Expectantly, BMI did not change progressively between the subgroups.

Figure 4

Classification of NAFL subgroups. A: Classification of NAFL subgroups against HO and NASH groups by ratios of SM 43:1;2 to SM 38:2;2 and SM41:3;2 to SM 38:2;2 whose trend lines are shown in Fig. 3C. B: Violin plots of log2-scaled abundances of CE 18:3, DG 34:1, and PE O- 38:7 in the three NAFL subgroups. Color coding of in panels A and B is the same as in Fig. 3C. C: Distribution of NAFL patients in bicluster components indicating the representative number of patients in each subgroup.

Next, we subjected the patients subgroups to random-forest classification according to the following scenarios: (a) h-NAFL versus HO, (b) hs-NAFL versus HO, (c) s-NAFL versus. HO, (d) s-NAFL versus NASH and (e) hs-NAFL versus NASH, (f) h-NAFL versus NASH on the full lipidome (supplemental Fig. S20). As anticipated, we observed that the classifications performance increased from (a) to (c) and from (d) to (f) indicating that the lipidome of NAFL subgroups consistently changed with the disease progression. The change of SM ratios in subgroups of NAFL and also in four major groups of patients visualized this trend (Fig. 3C).

In these classifications DG 34:1, but also CE 18:3 and CE 18:2 were the most discriminating markers (supplemental Fig. S21). We next tested, if h-NAFL, hs-NAFL, and s-NAFL patients, selected solely by comparing their lipidomes, are also distinguishable histologically? If so, what marker lipids could reflect the transition from h-NAFL to s-NAFL and further toward NASH? The examination of histological indices indicated gradual increase in both inflammation and steatosis (supplemental Table S6). However, none of them reached typical NASH values, apart from % of accumulated fat. In their lipidomes, the abundance of CE 18:3 and DG 34:1 gradually increased from h-NAFL to s-NAFL (Fig. 4B).

We also classified lipidomes of NAFL patients from the steatotic bicluster component against all remaining NAFL patients and identified PE O- 38:7 as the most discriminating marker (supplemental Fig. S22). Examination by HCD FT MS/MS identified it as plasmalogen PE O- 16:1/22:6 (supplemental Fig. S23 and supplemental Materials and Methods, MS2 validation) whose abundance dropped already in hs-NAFL subgroup down to s-NAFL level. Another three polyunsaturated glycerophospholipids: LPE 22:6, PC 38:6, and PC 40:7 showed similar classification power (supplemental Fig. S20). Their sum formula suggests that, similar to PE O- 38:7 and LPE 22:6, they also comprise C22:6 fatty acid moiety. We argue that stepwise decrease in the abundance of C22:6 - containing glycerophospholipids might indicate the onset of transition of h-NAFL toward s-NAFL, despite that there was only a marginal difference between histological indices of h-NAFL and hs-NAFL subgroups (supplemental Table S6). Interestingly, the classifications based on the SM ratios showed similar specificity (supplemental Fig. S17).

Taken together, our analysis distinguished nonsteatotic (NC and HO) and steatotic (NAFL and NASH) cohorts by the accumulation of storage lipid markers TG 50:1, DG 34:1, CE 18:3 and, independently, by the ratios of sphingomyelins SM 43:3;2 and SM 41:3;2 SM to the reference SM 38:2;2. Furthermore, the SM ratios together with coherently changing C22:6 - containing glycerophospholipids (most specifically PE O- 38:7, but also LPE 22:6, PC 38:6 and PC 40:7) revealed the compositional heterogeneity of the NAFL lipidome and how it changes, once the progressing disease gradually remodels HO-like lipidome toward NASH (Fig. 5).

Figure 5

Schematic trends of abundances of lipidomics and histological markers during NAFLD. The magnitude of change (in %) of lipid (Figs. 3 and 4) and histological (supplemental Table S4) markers scaled to the difference between NC and NASH. Cartoon images of the liver indicate progressive accumulation of fat eventually leading to inflammation and fibrosis.

Discussion

Excessive accumulation of TG is the histological hallmark of NAFLD. While this alone links NAFLD to altered lipid metabolism, little was known if and how NAFLD globally alters the liver lipidome, particularly its membrane and signaling complements.

Our study confirmed the accumulation of neutral lipids along with NAFLD progression. Interestingly, the molecular composition of TG and DG species did not change and, altogether, was similar to white adipose tissue (WAT). Accordingly, biclustering analysis identified monounsaturated TAG 50:1 and DAG 34:1, but also unsaturated CE 18:3 as most specific markers of NAFLD progression whose abundance steadily increases from HO to NASH (Figure 3, Figure 4, Figure 5). The impact of five risk genes implicated in NAFLD development corroborated the molecular specificity of PNPLA3 and MBOAT7, but did not alter the liver lipidome globally.

Interestingly, the levels of glycerophospholipids and of sphingolipids, e.g., Cer and (with a few notable exceptions) SM were not perturbed significantly (Fig. 1). Also, no apparent association with altered SM/Cer metabolism emerged from the transcriptomics analysis (69). Although liver inflammation and apoptosis in NASH are supposed to increase the levels of Cer also by cleaving SM (9, 13, 32, 70, 71), we did not observe it in the biopsies. However, if SM directly contributed to or are associated with steatosis and its transformation to NASH was an open question (72). Specific bidirectional changes of the four SM species supported the robust same lipid class disease state classification. Conveniently, ratios of SM markers could be computed from raw intensities of corresponding molecular peaks: they did not depend on the biopsy size and did not require the full lipidome quantification. To the best of our knowledge, the prospective marker lipids SM 43:3;2, SM 41:3;2, and SM 42:4;2 were not previously spotted in the pathophysiological context of metabolic syndrome. Therefore, it seems promising to use targeted quantification to follow their levels in plasma as it might yield a convenient marker for independent and noninvasive stratification of NAFLD.

Sphingolipids with n24:2 moities are produced from linoleic (C18:2) fatty acid (73). CE 18:2 is the most abundant cholesteryl ester in both the liver and plasma, and it is also enriched in NAFLD patients. It is therefore conceivable that enhanced biosynthesis of SM comprising n24:2 reflects increasing availability of free linoleic acid, while the abundance of SM with saturated or monounsaturated fatty acid moieties tends to decrease (Fig. 3B). However, relative changes of the abundance of SM species between steatotic and nonsteatotic conditions are not proportional to their absolute (molar) content and are also sphingosine-backbone dependent. We noticed that two out of three marker SM comprise sphingosine backbone having odd number of carbon atoms. Therefore, it would be interesting to assess the contribution of microbiome lipids (reviewed in (74)) or de novo synthesis from branched amino acids (75).

One of the most intriguing findings of this work was the lipidomics evidence of the NAFL cohort heterogeneity (Figure 3, Figure 4, Figure 5) that was not apparent from the histological examination (Table 1). In the NAFL subgroups that we conveniently termed as hs- and s-NAFL the levels of plasmalogen PE O- 38:7, but also PC 38:6, PC 40:7 and LPI 22:6 (all comprising C22:6 moiety) dropped down to the level typical for NASH patients. In the same NAFL subgroups the dynamics of histopathological indices (supplemental Table S6) was concordant with progressive steatosis with no hallmarks of fibrosis and inflammation. We therefore speculate that altered SM ratios together with decreased levels of C22:6-containing glycerophospholipids might indicate a turning point in NAFLD pathogenesis, where the transition from NAFL to NASH becomes irreversible.

While at this stage we cannot offer a suitable mechanistic explanation, we hypothesize that depletion of C22:6 containing GPL might affect the levels of PPARa (76) and SREBP1c (77) and, eventually, promote steatosis and inflammation. Also, we cannot rule out that the decreased abundance of these lipids hints at enhanced oxidative stress, despite that shotgun profiling revealed no notable accumulation of oxidized TG and glycerophospholipids.

Irrespective of NAFLD stage, the organism strived to maintain the compositional identity of the hepatocyte membrane lipidome and changes were limited in both molecular scope and magnitude. The unexpected link between species-specific SM metabolism and NAFLD progression, together with NAFL-specific change in the levels of C22:6 - containing glycerophospholipids could contribute to both patients stratification and mechanistic understanding of the lipidome regulation during the disease.

Last but not least, this work created an open and transparent lipidomic resource that could be expanded with or cross-validated by further independent studies reporting molar abundances of liver lipids, irrespective of their research objectives and employed analytical methods. Studies design and reported evidence should enable direct comparison of lipidomics data, rather than their context-dependent interpretations and trends. Eventually, this may help to establish reference values of molar lipid abundances and use them for molecular diagnostics of a broad spectrum of liver pathologies. It will also complement ongoing efforts to create harmonized lipidomic resources for liquid and solid biopsies sharing similar principles of collecting and organizing the data (41, 44, 55, 60, 78).

Data availability

All data concerned with this study are presented within the manuscript, supplemental data, and supplemental datasets.

Supplemental data

This article contains supplemental data (13, 17, 39, 42, 49, 50, 51, 64, 66, 67, 79, 80, 81, 82, 83, 84, 85, 86).

Conflict of interest

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