Lipidomic profiling identifies signatures of metabolic risk.

BACKGROUND: Metabolic syndrome (MetS), the clustering of metabolic risk factors, is associated with cardiovascular disease risk. We sought to determine if dysregulation of the lipidome may contribute to metabolic risk factors.
METHODS: We measured 154 circulating lipid species in 658 participants from the Framingham Heart Study (FHS) using liquid chromatography-tandem mass spectrometry and tested for associations with obesity, dysglycemia, and dyslipidemia. Independent external validation was sought in three independent cohorts. Follow-up data from the FHS were used to test for lipid metabolites associated with longitudinal changes in metabolic risk factors.
RESULTS: Thirty-nine lipids were associated with obesity and eight with dysglycemia in the FHS. Of 32 lipids that were available for replication for obesity and six for dyslipidemia, 28 (88%) replicated for obesity and five (83%) for dysglycemia. Four lipids were associated with longitudinal changes in body mass index and four were associated with changes in fasting blood glucose in the FHS.
CONCLUSIONS: We identified and replicated several novel lipid biomarkers of key metabolic traits. The lipid moieties identified in this study are involved in biological pathways of metabolic risk and can be explored for prognostic and therapeutic utility.

Research in context

Evidence before this study

Prior metabolomic association studies linking lipids with metabolic risk factors have been limited by the lipid metabolites profiled, lack of independent external replication, and absence of longitudinal association analyses.

Added value of this study

Our global lipidomics approach, along with external replication in three independent cohorts, provides robust cross-sectional associations of several lipid moieties with obesity and dysglycemia. We further strengthen these results by identifying lipids associated with longitudinal changes in metabolic risk factors.

Implications of all the available evidence

Our results, taken together with prior evidence, highlight lysoglycerophospholipids and dihydrosphingolipids as markers of obesity and dysglycemia that should be explored for prognostic and therapeutic utility.

Introduction

Cardiovascular disease (CVD), the leading cause of mortality in the United States and worldwide, accounts for roughly 30% of deaths [1] and is closely linked to the burden of metabolic risk factors. Metabolic syndrome (MetS) [2], clinically defined as the simultaneous presence of several metabolic risk factors [3], affects over 20% of American adults and is associated with a doubling in CVD risk [4]. While MetS is treated using lifestyle and drug-based interventions, its pathogenesis is not fully understood. Identifying key mechanistic drivers of metabolic risk factors and the progression from subclinical to overt disease may allow for the development of more targeted and effective interventions. The complexities of the associated metabolic risk factors necessitate global approaches, such as lipidomic profiling, which has advanced in the past decade due to mass spectrometry [5]. Considerable CVD research has shifted towards interrogating the circulating lipidome to further understand the cardiovascular consequences of dysregulated lipid metabolism [6].

Here, we report the results of plasma lipidomic profiling via liquid chromatography-tandem mass spectrometry (LC-MS/MS) in Framingham Heart Study (FHS) participants [7,8] in whom cross-sectional and longitudinal phenotypic characterization are available. Independent external replication of all cross-sectional associations from discovery was sought in the Erasmus Rucphen Family (ERF) study [9], San Antonio Family Heart Study (SAFHS) [10], and Progression of Early Subclinical Atherosclerosis (PESA) study [11]. We posited that the characteristic metabolic traits that define MetS have distinct and shared lipidomic signatures. To that end, we sought to identify lipid metabolomic biomarkers associated with metabolic risk factors to further understand the biology underlying MetS and to highlight high-value targets for the treatment and prevention of cardiometabolic disease.

Methods

Study samples

The discovery study sample included participants from the Offspring (n = 5124) and Third Generation (n = 4095) FHS cohorts (6,7). FHS clinical evaluations, including physical examinations and laboratory tests, were conducted every four to eight years. We sampled FHS Offspring cohort participants evaluated at Examination 8 (2005–2008; n = 3021) and Third Generation cohort participants evaluated at Examination 1 (2002–2005; n = 4095). Exclusion criteria included the following: off-site examination (n = 107), no blood specimen (n = 126), lack of unrestricted consent (n = 335), prevalent CVD at baseline (n = 216), pharmacological treatment for diabetes (n = 257) or dyslipidemia (n = 1201), missing BMI data (n = 4), cigarette smoking at the time of examination (n = 772), and age <50 or >79 years in the Offspring cohort (n = 123) or age <25 or >59 years in the Third Generation cohort (n = 183). The final study sample was selected from 3792 eligible participants.

Independent external replication was conducted in the ERF study, SAFHS, and PESA study. The ERF study is a prospective family-based study consisting of 3465 individuals aged 17–92 years from the southwest region of the Netherlands recruited based on kinship [9,12]. Samples without measures of metabolites or fasting glucose, or available genotype data for family kinship adjustment were excluded. The final study sample consisted of 2564 participants without diabetes. All baseline demographic data was collected between 2002 and 2006, and venous blood samples were collected after at least eight hours of fasting. A more detailed description of the ERF study and baseline characteristics of the study sample have been reported previously [9]

The SAFHS is a population-based, epidemiological study of type II diabetes (T2D) mellitus and CVD in Mexican Americans and non-Hispanic whites [10,13,14]. A total of 1431 individuals comprising 42 extended families who had given informed consent underwent blood sample collection. Of these, plasma samples from 1076 participants had complete data and were used for lipidomic profiling via liquid chromatography electrospray ionization tandem mass spectrometry (LC ESI-MS/MS) [13].

The PESA study is an observational cohort study of subclinical atherosclerosis in 4184 middle-aged (40–54 years) employees at Banco Santander in Madrid, Spain [11]. Recruitment began in 2010, and the baseline examination, including assessment for CVD risk factors and blood sampling, was completed in 2013. PESA participants with prevalent CVD at baseline (n = 34), active treatment for cancer (n = 2), history of transplants requiring active immunotherapy (n = 2), morbid obesity (BMI ≥40 kg/m2) (n = 22), chronic kidney disease (glomerular filtration rate ˂60 mL/min per square meter) (n = 1), an implanted metal device that could produce artifacts in computed tomography (n = 2), and any limitation or disability resulting in failure to comply with the study protocol (n = 1) were excluded. The final PESA sample size included 4120 participants.

Clinical measures

For the cross-sectional analyses, we evaluated three categorically-defined metabolic risk factors: obesity (BMI ≥30 kg/m2), dysglycemia (fasting blood glucose ≥100 mg/dL [adjusted for diabetes medications as described below]), and dyslipidemia (TAG ≥ 150 mg/dL [adjusted for lipid-lowering treatment as described below] or HDL cholesterol [HDL-C] <40 mg/dL in men or <50 mg/dL in women). All definitions were obtained from the National Cholesterol Education Program's Adult Treatment Panel III report [19]. For the longitudinal analyses, serial changes (from the baseline examination to the follow-up examination) in the continuous measures of the above-mentioned metabolic risk factors were used as the outcomes. The continuous measures for obesity, dysglycemia, and dyslipidemia were BMI [kg/m2], fasting blood glucose [mg/dL], and TAG [mg/dL] and HDL-C [mg/dL], respectively.

Study design

A 2 × 2 × 2 factorial study design was used in the FHS to group samples according to the presence or absence of each of the three central metabolic risk factors. Based on the clinical cut-points described above, eight (2 × 2 × 2 = 8) mutually-exclusive groups were generated from the three dichotomized outcomes. 80–84 individuals were selected for each group while controlling for sex and cohort across all groups. The longitudinal study sample consisted of all FHS participants in the cross-sectional sample who attended the subsequent clinic examination and had complete follow-up data for the traits of interest (n = 554). For those on lipid-lowering (n = 133) or diabetes medications (n = 30) at follow-up, observed values were also adjusted for treatment effects (see Statistical Methods for details). Fig. 1 provides a flow chart so the experimental design can be visualized.

Fig. 1

Flow chart of experimental design.

Lipidomic profiling

Lipidomic profiling for replication was carried out via three different platforms – 1) liquid chromatography-mass spectrometry (LC-MS), 2) AbsoluteIDQ p150 Kit of Biocrates Life Sciences AG, and 3) electrospray-ionization mass spectrometry (ESI-MS) – using methods previously described [15], [16], [17], [18]

A detailed description of the LC-MS based lipidomic profiling methods used in the FHS is provided in the Supplementary Text and Table 2a contains the metabolite transitions used in the MRM mass spectrometry.

Table 2

Cross-sectional lipid associations with obesity in FHS.

Lipid species	Class	Estimated β*	SE	P-value†	Replication#
LPC 17:0	LPC	−0.618	0.072	1.05E-16	✓
LPC 18:2	LPC	−0.524	0.069	1.67E-13	✓
LPC 18:0p	LPC	−0.503	0.070	2.28E-12	✓
LPC 18:1	LPC	−0.510	0.072	3.89E-12	✓
LPC 18:0e	LPC	−0.503	0.071	3.95E-12	✓
PC 16:2e	PC	−0.448	0.067	4.52E-11
SM 36:0	SM	0.493	0.075	7.81E-11
LPC 22:5	LPC	−0.462	0.073	5.08E-10
PC 40:8	PC	−0.460	0.074	8.80E-10
LPC 17:1	LPC	−0.447	0.073	1.97E-09	✓
LPC 15:0	LPC	−0.439	0.076	1.06E-08	✓
LPC 18:0	LPC	−0.423	0.074	1.34E-08	✓
LPC 22:6	LPC	−0.427	0.075	2.26E-08	✓
PC 34:0	PC	−0.407	0.074	6.45E-08
LPC 16:0	LPC	−0.387	0.072	1.31E-07	✓
LPC 20:4	LPC	−0.365	0.071	3.30E-07	✓
LPC 18:1e	LPC	−0.373	0.073	4.35E-07	✓
TAG 50:1	TAG	0.329	0.065	6.44E-07	✓
TAG 50:2	TAG	0.323	0.065	7.76E-07	✓
LPC 18:3	LPC	−0.369	0.076	1.40E-06
PC 38:3	PC	0.362	0.075	1.66E-06	✓
SM (d18:0/24:0)	SM	0.347	0.075	4.86E-06
LPE 16:0	LPE	−0.352	0.078	7.17E-06	✓
TAG 48:0	TAG	0.310	0.069	8.65E-06	✓
LPE 18:1	LPE	−0.343	0.077	1.04E-05	✓
LPC 16:0e	LPC	−0.317	0.072	1.16E-05	✓
LPC 16:1e	LPC	−0.310	0.072	1.86E-05	✓
DAG 32:0	DAG	0.280	0.067	3.04E-05
LPC 18:2e	LPC	−0.299	0.074	5.24E-05	✓
TAG 48:1	TAG	0.276	0.069	6.40E-05	✓
PC 40:7	PC	−0.302	0.075	6.85E-05	✓
PC 36:3e	PC	−0.283	0.072	1.01E-04	✓
LPC 20:5	LPC	−0.296	0.076	1.04E-04
TAG 52:1	TAG	0.251	0.066	1.51E-04	✓
LPE 18:2	LPE	−0.276	0.078	4.26E-04	✓
PC 34:1e	PC	−0.249	0.071	4.58E-04
PC 36:2	PC	−0.270	0.077	4.85E-04	✓
PC 36:0	PC	−0.270	0.077	5.13E-04	✓
DAG 34:1	DAG	0.221	0.064	5.48E-04

Estimated β coefficients represent the mean differences in standardized lipid measures between participants with and without obesity.

The p-value threshold for significance was determined by the Bonferroni method and based on the number of principal components. In FHS, the p-value threshold was p < 5.7 × 10−4 [0.05/(22×4)].

The p-value threshold for replication in the ERF study, PESA cohort, and SAFHS of the lipid cross-sectional associations with obesity was calculated using Bonferroni correction to be 1.56 × 10−3 (0.05/32).

✓ denotes lipid species that replicated in the ERF study, PESA cohort, or SAFHS.

Abbreviations: DAG=diacylglycerol; LPC=lysophosphatidylcholine; LPE=lysophosphatidylethanolamine; PC=phosphatidylcholine; SM=sphingomyelin; TAG=triacylglycerol.

Statistical methods

Statistical analyses were performed with SAS software version 9.2 (SAS Institute, Inc., Cary, NC, USA). All lipid measurements were rank-normalized via the Blom method [20]. For the FHS cross-sectional analyses, we used generalized linear mixed models (SAS procedure: GLIMMIX) to identify lipid species that (a) differed across all eight metabolic risk factor groups and (b) were associated with individual risk factors. All analyses were adjusted for the following baseline covariates: age, sex, cohort, and age-cohort interaction. To identify lipid-level associations across all eight groups, each lipid marker was treated as the outcome variable and the eight-level risk group as the main predictor. That is, our models were lipid level ~ cohort age sex cohort*age risk-factor-classes (* levels).

To identify independent associations between lipid levels and metabolic risk factors, each lipid marker was considered as the outcome variable and each risk factor as the main predictor with the other two risk factors as covariates. That is, our models were lipid level ~ cohort age sex cohort*age obesity dysclycemia dyslipidemia, where cohort, sex, obesity, dysclycemia and dyslipidemia are binary variables.

Secondary analyses of the dyslipidemia axis were conducted using the same models to explore independent associations between lipid species and high versus low TAG and HDL-C. In the PESA samples, generalized linear mixed models were used to identify lipid associations across all eight risk factor groups and with obesity, dysglycemia, and dyslipidemia. All replication models were adjusted for sex and baseline age. In addition, for individuals in the PESA cohort with self-reported use of lipid-lowering or diabetes medications, values were adjusted for treatment effects as follows. HDL-C was lowered by 10% and TAG was raised by 20% for individuals on lipid-lowering treatment, and glucose was raised by 10% or to a maximum of 126 mg/dL for those on diabetes medications [21]. In the ERF, pairwise partial correlation coefficients were generated for each metabolite-trait association, and all analyses were adjusted for age, sex, lipid-lowering medications, and family relationship using polygenic in the GenABEL R package [22]. Bonferroni-corrected p-values for obesity and dysglycemia were individually calculated as 0.05/n where n was equal to the number of metabolites measured in the ERF with available data for the corresponding trait.

For the longitudinal analyses, single marker associations were evaluated using general linear models (SAS procedure: GLM) with each lipid marker as the predictor and longitudinal changes in each metabolic trait (BMI, glucose, TAG, and HDL-C) as the outcome. All models were adjusted for age, sex, cohort, and baseline values of BMI, glucose, TAG, total cholesterol, and HDL-C. For individuals with self-reported use of lipid-lowering or diabetes medications at follow-up, values were adjusted as described above. We evaluated multimarker associations with longitudinal changes using stepwise regression with forward selection. Lipid markers with false discovery rate (FDR) < 0.25 in the single-marker analysis were considered as initial candidates for each trait. A dynamic p-value threshold corresponding to FDR < 0.25 was used for selection at each step [at step i, p < (0.25*i/154)]. For each trait, an aggregate R2 was calculated from the final model after forward selection.

To account for intercorrelations among lipid markers, principal component analyses were performed for the discovery (n = 22 principal components) and PESA replication (n = 13 principal components) sets. The number of principal components with eigenvalues >1 was deemed the number of “independent” tests, and Bonferroni correction was applied to account for multiple testing. For the cross-sectional analyses, the experiment-wide p-value threshold was determined as 0.05/n where n was equal to the number of principal components times four primary analyses (three risk factors plus the eight-group comparison) [23]. This method yielded a threshold of p < 5.7 × 10−4 for the FHS [0.05/(4 × 22)]. The same methodology was used to determine the p-value threshold (p < 5.7 × 10−4) for the longitudinal study, considering the four primary analyses (one per continuous metabolic trait).

Ethics statement

Informed consent was given by all participants at time of exam cycle from the FHS. The clinical study protocol was approved by the Boston University Medical Campus Institutional Review Board.

Results

Study sample characteristics

Table 1 lists the clinical characteristics of the 658 FHS participants according to the eight prespecified combinatorial groups. By design, mean age and percentage of women were consistent across all eight groups in the FHS, while by design, BMI, glucose, and lipid levels (HDL-C and TAG) varied. Supplementary Table 1 shows the clinical characteristics of the FHS longitudinal sample (n = 554) at the baseline and follow-up examinations. Mean changes, adjusted for lipid-lowering or hyperglycemia medication-use, were determined from baseline to follow-up for the four metabolic traits: 0.6 ± 2.6 kg/m2 for BMI, 0.5 ± 20.4 mg/dL for fasting glucose, 3.6 ± 10.7 mg/dL for HDL-C, and −0.5 ± 84.6 mg/dL for TAG.

Table 1

Study design and clinical characteristics of FHS samples.

Clinical characteristics	Group 1	Group 2	Group 3	Group 4	Group 5	Group 6	Group 7	Group 8
Obesity*	–	–	–	–	+	+	+	+
Dyslipidemia*	–	–	+	+	–	–	+	+
Dysglycemia*	–	+	–	+	–	+	–	+
Sample size (n)	84	83	84	84	82	84	74	83
Age (years)	52 ± 14	54 ± 13	52 ± 14	54 ± 15	53 ± 12	54 ± 11	50 ± 14	55 ± 13
Women (%)	50	49	50	50	50	50	57	49
BMI (kg/m2)	24.5 ± 2.8	25.0 ± 2.7	25.9 ± 2.7	26.6 ± 2.3	32.4 ± 2.3	35.9 ± 5.1	35.4 ± 4.8	34.9 ± 4.9
HDL-C (mg/dL)	65 ± 18	62 ± 15	45 ± 11	44 ± 14	61 ± 13	58 ± 13	42 ± 10	42 ± 10
TAG (mg/dL)	81 ± 26	84 ± 26	167 ± 115	175 ± 91	92 ± 28	95 ± 29	165 ± 77	202 ± 148
Glucose (mg/dL)	92 ± 5	107 ± 15	92 ± 5	107 ± 7	93 ± 4	109 ± 7	93 ± 5	122 ± 37

(+) and (-) represent the presence and absence, respectively, of obesity, dysglycemia, and dyslipidemia in each of the eight metabolic risk factor groups. Obesity was defined as BMI ≥30 kg/m2. Dysglycemia was defined as fasting blood glucose ≥100 mg/dL. Dyslipidemia was defined as either TAG ≥150 mg/dL or HDL-C <40 mg/dL in men or <50 mg/dL in women.

Analytical set of potential lipidomic biomarkers

154 distinct metabolites representing ten classes and subclasses were measured in the FHS: four ceramides (Cers), 11 cholesteryl esters (CEs), eight diacylglycerols (DAGs), 21 lysophosphatidylcholines (LPCs), four lysophosphatidylethanolamines (LPEs), 28 phosphatidylcholines (PCs), 11 phosphatidylethanolamines (PEs), one phosphatidylserine (PS), 26 sphingomyelins (SMs), and 40 triacylglycerols (TAGs). A correlation matrix for all 154 lipid species is provided in Supplementary Figure 1. Many metabolite species were positively correlated, with the exception of some LPCs and TAGs, which were negatively correlated. Fig. 1 provides a flow chart of the analytical set and how each marker are used for analysis. .

FHS cross-sectional associations with metabolic risk factors

Among the aforementioned lipid species, several were cross-sectionally associated with obesity (n = 39; Table 2), dysglycemia (n = 8; Table 3), and dyslipidemia (n = 89; Supplementary Table 4) in the FHS. Many of the lipid markers for obesity were LPCs (n = 18; Table 2). LPCs, PCs, and LPEs were inversely associated, while SMs, TAGs, and DAGs were directly correlated with obesity (Table 2). Eight lipids were associated with dysglycemia; LPCs (n = 5, Table 3) were inversely associated, while SMs and Cer (n = 3, Table 3) were directly correlated with dysglycemia. The SM and Cer species are dihydrosphingolipids (dhSLs), which have backbones comprised of sphinganine (d18:0) rather than sphingosine (d18:1) as is the case for sphingolipids (SLs). Supplementary Table 4 lists all 89 lipid species associated with dyslipidemia (direct = 63, inverse = 26) along with their corresponding replication results. We further analyzed the metabolites by class; results are shown in Supplementary Figures 2, 3, 4, and 5 for obesity, dysglycemia, low vs. high HDL-C, and high vs. low TAG, respectively. Despite directional consistency, the majority of the associations between metabolite classes representing several individual lipid species and metabolic risk factors were not statistically significant, likely due to the broad diversity of individual metabolite species represented.

Table 3

Cross-sectional lipid associations with dysglycemia in FHS.

Lipid Species	Lipid Class	Estimated β*	SE	P-value†	Replication#
Cer (d18:0/24:0)	Cer	0.280	0.070	6.59E-05
LPC 18:1e	LPC	−0.279	0.073	1.50E-04	✓
LPC 17:0	LPC	−0.271	0.073	1.98E-04	✓
SM 36:0	SM	0.275	0.075	2.43E-04
LPC 18:0e	LPC	−0.260	0.071	2.80E-04	✓
LPC 16:1e	LPC	−0.254	0.072	4.48E-04	✓
SM (d18:0/24:0)	SM	0.265	0.075	4.79E-04
LPC 18:0p	LPC	−0.244	0.070	5.48E-04	✓

Estimated β coefficients represent the mean differences in standardized lipid measures between participants with and without dysglycemia.

The p-value threshold for significance in FHS (p < 5.7 × 10−4) was determined by the Bonferroni method [0.05/(22×4)].

The p-value threshold for replication in the ERF study, PESA cohort, and SAFHS of the lipid cross-sectional associations with dysglycemia was calculated using Bonferroni correction to be 8.33 × 10−3 (0.05/6).

✓ denotes lipid species that replicated in the ERF, SAFHS, or PESA.

Abbreviations: Cer = ceramide; LPC = lysophosphatidylcholine; SM = sphingomyelin.

External replication

Independent external replication of all FHS cross-sectional metabolite associations with obesity and dysglycemia was sought in the ERF study, SAFHS, and PESA. Of the 39 metabolites that were cross-sectionally associated with obesity in the FHS, a total of 32 had available data for replication. Of these 32 metabolites, 28 were significantly associated with obesity at p < 1.56E-03 (0.05/32) (Table 2, Supplementary Table 2). Of the eight metabolites cross-sectionally associated with dysglycemia in the FHS, six had available data for replication. Of these six metabolites, five were associated with dysglycemia at p < 8.3E-03 (0.05/6) (Table 3, Supplementary Table 3).

Associations with longitudinal changes in metabolic traits

We identified several lipid species associated with longitudinal changes in BMI (n = 4), glucose (n = 4), and HDL-C (n = 1) (Table 4). LPCs 16:0, 18:0, 16:0e, and 22:5 were inversely associated with change in BMI [follow-up minus baseline]. These four LPCs were also inversely associated with BMI in the cross-sectional analyses. LPC 18:1 and PC 36:2 were inversely associated with change in glucose, and SM 36:0 and SM (d18:0/24:0) were directly associated with change in glucose. Both SMs were also directly associated with dysglycemia in the cross-sectional analyses. SM (d18:2/24:1), which was inversely associated with dyslipidemia in cross-sectional analysis, was directly associated with longitudinal change in HDL-C.

Table 4

Single marker associations with longitudinal changes in metabolic risk factors.

Change in metabolic trait	Lipid species	Lipid class	Estimated β*	SE	P-value†	FDR	Partial R2
BMI	LPC 16:0	LPC	−0.57	0.14	6.7E-05	0.01	0.029
	LPC 18:0	LPC	−0.50	0.13	1.7E-04	0.01	0.026
	LPC 16:0e	LPC	−0.44	0.12	4.2E-04	0.02	0.023
	LPC 22:5	LPC	−0.45	0.13	4.3E-04	0.02	0.023
Glucose	SM 36:0	SM	3.71	0.91	5.0E-05	0.01	0.030
	LPC 18:1	LPC	−3.37	0.93	2.9E-04	0.02	0.024
	SM (d18:0/24:0)	SM	3.34	0.92	3.2E-04	0.02	0.024
	PC 36:2	PC	−2.97	0.85	5.2E-04	0.02	0.022
HDL-C	SM (d18:2/24:1)	SM	2.61	0.55	2.8E-06	0.00	0.040

Estimated β coefficients represent the mean differences in longitudinal change (native units) in the stated metabolic trait per one standard deviation difference in lipid measures.

The p-value threshold for significance (p < 5.7 × 10−4) was determined by the Bonferroni method [0.05/(22×4)].

Abbreviations: LPC = lysophosphatidylcholine; PC = phosphatidylcholine; SM = sphingomyelin.

Lipids also associated in cross-sectional analysis with the corresponding metabolic risk factor are shown in bold.

In multimarker analyses, we identified multimarker panels associated with longitudinal changes in metabolic traits (Table 5). The BMI multimarker consisted of LPC 16:0, which explained 2.9% of the variation in longitudinal change in BMI (baseline BMI explained 2.7%). SM 36:0 and DAG 36:2 were retained in the final glucose model, together accounting for 6.1% of the variation in change in glucose (baseline glucose explained 22.1%). The HDL-C model contained SM (d18:2/24:1), LPC 20:3, and PS 36:0. This multimarker panel explained 8.4% of the variation in change in HDL-C (baseline HDL-C explained 3.2%). LPC 18:3, SM (d18:1/16:0), and Cer (d18:1/24:1) were retained in the final TAG model. These lipids accounted for 8% of the variation in change in TAG (baseline TAG explained 20.7%).

Table 5

Multimarker associations with longitudinal changes in metabolic risk factors.

Change in metabolic trait	Parameter	Estimated β*	SE	P-value†	R[2]
BMI	Baseline BMI	−0.08	0.02	1E-04	0.027
	Baseline Covariates‡				0.051
	LPC 16:0	−0.56	0.14	<1.0E-04	0.029
	Multi-Metabolite Panel				0.029
Glucose	Baseline Glucose	−0.53	0.04	<1.0E-04	0.221
	Baseline Covariates‡				0.289
	SM 36:0	4.05	0.90	<1.0E-04	0.036
	DAG 36:2	−4.58	1.23	2.0E-04	0.025
	Multi-Metabolite Panel				0.061
HDL-C	Baseline HDL-C	−0.15	0.04	<1.0E-04	0.032
	Baseline Covariates‡				0.097
	SM (d18:2/24:1)	2.67	0.57	<1.0E-04	0.042
	LPC 20:3	1.60	0.50	1.6E-03	0.020
	PS 36:0	−1.60	0.48	9.0E-04	0.022
	Multi-Metabolite Panel				0.084
TAG	Baseline TAG	−0.57	0.05	<1.0E-04	0.207
	Baseline Covariates‡				0.260
	LPC 18:3	−12.75	3.61	4.0E-04	0.023
	SM (d18:1/16:0)	−21.18	5.05	<1.0E-04	0.032
	Cer (d18:1/24:1)	15.91	4.24	2.0E-04	0.026
	Multi-Metabolite Panel				0.080

Estimated β coefficients represent the mean differences in longitudinal change (native units) in the stated metabolic trait per one standard deviation difference in lipid measures.

A dynamic p-value threshold for entrance into the model was defined based on FDR < 0.25 [at step i, p < (0.25*i)/154].

Baseline covariates were as follows: age; sex; cohort; batch; and baseline values for total cholesterol, HDL-C, glucose, and BMI.

Abbreviations: LPC = lysophosphatidylcholine; DAG = diacylglycerol; PS = phosphatidylserine; SM = sphingomyelin; Cer = ceramide.

Lipids also significant in individual longitudinal associations are shown in italics.

Discussion

Associations between plasma lipid classes and metabolic risk factors have been reported [24], [25], [26], [27], [28], [29], but many of these studies were limited to single site studies with small numbers of lipid species or case-control analyses. We used liquid chromatography-tandem mass spectrometry to detect cross-sectional and longitudinal lipidomic signatures of key metabolic risk factors. Our cross-sectional lipidomic signatures of obesity and dysglycemia were highly replicable in external studies and highlight two lipid classes as biomarkers of metabolic risk: lysoglycerophospholipids (LGPLs) and sphingolipids (SLs).

Lysoglycerophospholipid markers of metabolic traits

Cross-sectional and longitudinal analyses of obesity/BMI and dysglycemia/fasting blood glucose identified significant associations with LGPLs, which are generated along with a free fatty acid by phospholipase-mediated hydrolysis of glycerophospholipids (GPLs; supplementary figure 6c) [30]. The predominant LGPL in our results was lysophosphatidylcholine (LPC). We identified inverse cross-sectional associations between all LPCs and obesity (n = 18; Table 2). Fifteen (83%) of these associations replicated in the ERF and SAFHS (Supplementary Tables 2 and 3). Four of the 18 LPCs cross-sectionally associated with obesity in the FHS were also associated with longitudinal changes in BMI (Table 4). Five of the LPCs associated with obesity cross-sectionally in the FHS were also inversely associated with dysglycemia (Supplementary Table 5), and all five of these associations replicated. Nine phosphatidylcholines (PCs) were associated with obesity (Table 2) – all showing inverse associations, with the exception of PC (38:3), which showed a positive association with obesity (β = 0.36, P = 1.7E-06; Table 2). This association replicated in the ERF, and prior work has shown PC (38:3) to be positively associated with CVD [31] and enriched in atherosclerotic plaques [32]. Of the nine PCs associated with obesity, seven were available for replication and five (71%), including PC (38:3), replicated (Supplementary Table 2). Lysophosphatidylethanolamine (LPE), another LGPL class detected by our lipid-profiling platform, was inversely associated with obesity in the FHS (Table 2); all three LPEs associated with obesity in the FHS replicated (Supplementary Table 2). Overall, LGPLs were inversely associated with metabolic risk factors.

Recent evidence, which has implicated LGPLs in essential signaling processes [33], highlights the potential utility of these lipid classes as novel biomarkers of cardiometabolic risk. LPCs are the most well-characterized and abundant LGPL in human plasma (125–143 µM in healthy individuals) [34], and circulating levels are determined by their generation via two pathways – 1) transfer of a fatty acid from PC to cholesterol by lecithin-cholesterol acyltransferase (LCAT; supplementary Figure 6d), and 2) phospholipase A2 (PLA2)-mediated hydrolysis of PC [35] (Supplementary Figure 6c, bottom) and clearance via transporters such as albumin and alpha-1 acid glycoprotein [34]. While LPCs have been found to be elevated in atherosclerotic plaques [32], and thus might be expected to positively correlate with metabolic risk factors, we found inverse associations between LPCs and obesity and dysglycemia (Tables 2, 3). This seemingly counterintuitive finding for LPCs, however, has been demonstrated for CVD [31,36] and perhaps reflects an increase in LPC catabolism or a decrease in its generation in such physiological states. Indeed, individuals with coronary artery disease have been shown to have lower LCAT activity [37], and three recent studies found reduced plasma LPC levels in obese patients [24,25,13]. Furthermore, higher levels of circulating LPC could reflect a lower rate of atherogenic lysophosphatidic acid production via autotaxin-mediated LPC hydrolysis, thus conferring a net protective effect on CVD [38]. Therefore, in conjunction with recent epidemiological evidence associating dairy consumption with changes in circulating LPC 17:0 [39], an abundance of genetic, epidemiologic, functional, and biochemical data implicating LPCs in CVD suggest that mechanistically understanding the enzyme kinetics and complexities of LPC metabolism may highlight promising therapeutic targets.

Sphingolipid markers for metabolic traits

Sphingolipids, comprised of Cers, SMs, and their reduced forms (dhCers and dhSMs), were positively associated with obesity (Table 2) and dysglycemia (Table 3). The two dhSMs positively associated with obesity – SM 36:0 and SM (d18:0/24:0) – were also positively associated with dysglycemia, along with Cer (d18:0/24:0). Cer (d18:0/24:0) shares the same sphingosine backbone and fatty acyl side chains as SM (d18:0/24:0) such that the two are inter-convertible via the enzymatic actions of sphingomyelin synthase (SMS; Cer→SM) and sphingomyelinase (SMase; SM→Cer). Of note, SM 36:0 and SM (d18:0/24:0) were also associated with longitudinal change in glucose (Table 4).

Sphingolipids, which are characterized by a sphingosine backbone, are synthesized de novo from the condensation of palmitoyl-CoA and serine by serine palmitoyltransferase [40], with the resulting dhCer being oxidizable by dihydroceramide desaturases (Des) to Cer. Sphingolipids are integral structural constituents of lipid membranes and are involved in cell signaling and apoptosis [41,42], with SMs being the most prevalent SLs in human plasma that can be derived from Cers via SMS, releasing a DAG in the process [41,43] (Supplementary Figure 6a). Sphingomyelin inhibits lipoprotein lipase (LPL)-catalyzed hydrolysis of TAG [44] and LCAT-catalyzed esterification of cholesterol [45] (Supplementary Figure 6b), and individuals with SMase deficiency (Niemann-Pick disease) have low plasma HDL-C levels [46] due to SM-mediated inhibition of cholesterol efflux. Sphingomyelin accumulation has been associated with T2DM [47] and observed in mouse models of insulin resistance [48,49] as well as in the skeletal muscle and serum of obese, insulin-resistant humans [27,50,51].

Prior studies have shown plasma SM levels to predict coronary artery disease risk, possibly due to their strong correlations with apolipoprotein B (apoB) and TAG levels [52]. Indeed, the SM content of apoB-containing or TAG-rich lipoproteins is higher than that in HDL [53], and hydrolysis of LDL-associated SM to Cer causes LDL particle aggregation and macrophage foam cell formation [54,55], which are two hallmarks of atherosclerosis. Recent studies have suggested that Cer provides additional prognostic utility for CVD beyond that from LDL-C [56] – perhaps because it acts as a proxy to standardize changes in LDL-C mass by changes in apoB concentration, which is the biological determinant of atherogenicity of a high LDL-C level [57]. Interestingly, a missense mutation (rs267738) in the CERS2 gene, which encodes the enzyme ceramide synthase that is responsible for biosynthesis of very long-chain Cer (C22/C24/C24:1), is also associated with hemoglobin A1c levels [58], and functional studies have shown that CerS2 haploinsufficiency in mice results in insulin resistance and non-alcoholic fatty liver disease [59], possibly due to a compensatory increase in atherogenic Cer-16. Of note, several dietary intervention studies have demonstrated that improved diet quality blunts the association between Cer and CVD [60], perhaps due to the dependence of Cer synthesis on dietary intake of saturated fats [61]. Nonetheless, the most direct evidence comes from mouse knockout studies in which genetic ablation of key enzymes in Cer synthesis mitigates atherosclerosis [62], insulin resistance [63], and T2DM [64] among other metabolic-related disease phenotypes. Taken together with our results, there is an abundance of evidence for Cers and SMs as markers of obesity and dysglycemia risk.

We acknowledge several limitations of our study. First, the effect sizes and metabolite species differed among the FHS and replication cohorts, despite directional consistency in the replicable metabolites. These discrepancies were in part due to differences in environmental exposures (i.e. smoking) [13,65] as well as discrepancies in lipidomic and analytical techniques across studies. Moreover, some lipids were unavailable for replication. Even considering these limitations, however, consistent relationships were observed, thus reinforcing the robustness of the replicated results. Second, the longitudinal changes in mean TAG levels (<1 mg/dL over six-years of follow up) were small and not of a magnitude to confer clinical consequences. Third, our lipid quantification methods did not allow for the separation of isomeric species and thus multiple isomers could have contributed to each lipid measurement. Fourth, the low replication rate for sphingolipids was partially due to the fact that many of the sphingolipids quantified in the FHS via liquid chromatography tandem mass spectrometry were not measured in the replication cohorts, which used different lipidomic platforms. Finally, larger-scale, multi-site epidemiological studies as well as functional experiments are needed to confirm and extend the findings of this investigation.

Conclusions

We identified numerous lipid species that were associated with metabolic risk factors cross-sectionally and with their longitudinal changes. Of note, we demonstrated robust and replicable associations of lysoglycerophospholipids and dihydrosphingolipids with obesity and dysglycemia. Taken together with biological evidence and prior interventional and functional studies, we hypothesize that the lipid species identified in this study may serve as key mechanistic drivers of metabolic dysregulation and can lay the groundwork for future studies to explore their potential therapeutic relevance.

CRediT authorship contribution statement

Xiaoyan Yin: Conceptualization, Writing - original draft. Christine M. Willinger: Conceptualization, Writing - original draft. Joshua Keefe: Data curation, Formal analysis, Writing - review & editing. Jun Liu: Project administration, Validation. Antonio Fernández-Ortiz: Project administration, Validation. Borja Ibáñez: Project administration, Validation. José Peñalvo: Project administration, Validation. Aram Adourian: Project administration, Validation. George Chen: Project administration, Validation. Dolores Corella: Project administration, Validation. Reinald Pamplona: Project administration, Validation. Manuel Portero-Otin: Project administration, Validation. Mariona Jove: Project administration, Validation. Paul Courchesne: Data curation, Project administration. Cornelia M. van Duijn: Project administration, Validation. Valentín Fuster: Project administration, Validation. José M. Ordovás: Data curation, Formal analysis. Ayşe Demirkan: Project administration, Validation. Martin G. Larson: Data curation, Formal analysis. Daniel Levy: Conceptualization, Supervision, Writing - review & editing.

Declaration of Competing Interest

Dr. Adourian contributed to this project while affiliated with BG Medicine, Inc., a biomarker discovery company. All other authors have reported that they have no relationships relevant to the contents of this paper to disclose.