In this retrospective study, a cohort of 106 adult patients (74 men, 32 women, median age 49, 95% CI 49–53 years) with chronic heart failure (i.e., at least one previous heart failure event, comprising hospitalization for HF and/or an urgent visit resulting in intravenous therapy for HF due to DCM) was examined. The present cohort is a sub-group of the population analyzed in our previous publication (28). The HF has been defined as a clinical syndrome characterized by fundamental symptoms (e.g., breathlessness, ankle swelling, and fatigue) and/or signs (e.g., elevated jugular venous pressure, pulmonary crackles, and peripheral edema) related to a structural and/or functional abnormality of the heart that results in elevated intracardiac pressures and/or inadequate cardiac output at rest and/or during exercise (34). A period of clinical stability of at least 6 months in optimal medical therapy (OMT) was required for enrollment. Patients were diagnosed between 1982 and 2011 (t0). In the period 2010–2012 (t1), patients were re-examined or examined for the first time (only for patients diagnosed between 2010 and 2012), blood serum samples were collected and analyzed via NMR. Then, survival status was evaluated at the last follow-up in 2019 (t2). In patients who died or were transplanted, the end of follow-up was considered either the time of death or heart transplantation. In the minority of patients lost to follow-up (i.e., not traceable by June 2019), the last clinical evaluation or the last telephone contact was considered. The experimental design is graphically illustrated in Figure 1.

The enrolled patients were classified as idiopathic DCM, defined by the presence of left ventricular (LV) or biventricular dilatation and systolic dysfunction in the absence of abnormal loading conditions (hypertension, valve disease) or coronary artery disease sufficient to cause global systolic impairment (1, 35). The patients with HF, judged to be secondary to ischemic heart disease, systemic hypertension, chemotherapy, alcoholic abuse, diabetes mellitus, cor pulmonalis, valve disease, or other cardiac or systemic diseases, were excluded, as well as patients for whom a coronary angiogram was not available. The patients were consecutively enrolled in the years 2010–2012 at the Careggi University Hospital Florence, Italy. They underwent half-yearly clinical assessments (median follow-up from enrollment 8 years). All study patients were evaluated and followed up by clinical history, physical examination, 12-lead ECG, standard chest radiograph, routine laboratory tests, M-mode, 2D, and Doppler echocardiography. For the entire study period, the patients were seen by the same cardiologists who assumed primary responsibility for their management.

This study was approved by the local Ethics Committee (Azienda Ospedaliero—Universitaria Careggi, Florence, Italy). Written informed consent was obtained from each participant at the time of blood sample collection. The study adheres to the principles of the Helsinki Declaration and its later amendments.

Each blood sample was collected through peripheral venous access in a 10-mL tube (BD P100, BD Diagnostics, Franklin Lakes, NJ). Subsequently, blood samples were centrifuged for 10 min at 4,000 rpm at the temperature of 4°C, then, the supernatant serum was aliquoted in sterile cryovials and stored at −80°C pending NMR analysis.

Serum samples were prepared for NMR experiments as described in our previous publications (28, 36). One-dimensional ¹H NMR spectra of each sample were acquired using a Bruker 600 MHz spectrometer (Bruker BioSpin) operating at 600.13 MHz proton Larmor frequency and equipped with a 5 mm CPTCI ¹H-¹³C-³¹P and ²H-decoupling cryoprobe, including a z-axis gradient coil, an automatic tuning-matching, and an automatic sample changer. A BTO 2000 thermocouple served for temperature stabilization at the level of ~0.1 K at the sample. Before measurement, samples were kept for at least 3 min inside the NMR probe-head, for temperature equilibration at 310 K. For each serum sample, a standard nuclear Overhauser effect spectroscopy pulse sequence NOESY 1Dpresat (noesygppr1d.comp; Bruker BioSpin)(37), using 64 scans, 98,304 data points, a spectral width of 18,028 Hz, an acquisition time of 2.7 s, a relaxation delay of 4 s and a mixing time of 0.01 s (total duration of the NMR experiment 7 min), was applied to obtain a spectrum in which both signals of low molecular weight metabolites and high molecular weight macromolecules (i.e., proteins, lipids, and lipoproteins) are detected.

Free induction decays were multiplied by an exponential function equivalent to 1.0 Hz line-broadening factor before applying Fourier transform. The transformed spectra were automatically corrected for phase and baseline distortions and calibrated to the anomeric glucose doublet at δ 5.24 ppm using TopSpin3.6.2 (Bruker Biospinsrl).

Each 1D spectrum was segmented into 0.02 ppm chemical shift bins in the range 0.2–10 ppm, and the corresponding spectral areas were integrated using the AssureNMR software (Bruker BioSpin). The regions of residual water signal (4.37–5.13 ppm) and the signals of ethanol (1.12–1.23 ppm and 3.53–3.73 ppm) were removed, and the dimension of the system was reduced to 438 bins. The probabilistic quotient normalization (38) was applied on the remaining bins prior to statistical analysis.

Data analyses were performed using the open-source software R. Multivariate analysis was performed on binned spectra (thus, on the whole spectra, considering both assigned and unassigned metabolites). Principal component analysis (PCA) was used as the first unsupervised analysis to visualize data. Metabolomics analysis was performed using a fingerprinting approach: the metabolomic fingerprint is a global, rapid evaluation of an NMR spectrum as a whole that considers all (assigned or unassigned) detectable metabolites present in that biological sample (39). Standard partial least square discriminant analysis (40) (PLS-DA) was applied to discriminate the metabolomic fingerprints of survivors and deceased patients using the first 7 PLS components, and the PLS-DA model was validated using a Leave-One-Out cross-validation scheme (LOOCV, R script developed in-house). Since the group size is unbalanced (survived 75.5%, deceased 24.5%), samples from 25 survivors and 25 deceased patients were randomly chosen from the full dataset and subjected to PLS-DA modeling. The resampling procedure was performed 100 times to account for variability in the sampling procedure, and each model was cross-validated each time. Samples were assigned to one of the two classes using the majority vote algorithm (R library “mclust”) on the results obtained by the 100 iterations. Sensitivity, specificity, and accuracy were calculated according to the standard definitions. Variable importance in projections (VIP) was calculated using an R script developed in-house, variables with a VIP score higher than 1 were considered important in the PLS-DA model.

The predictive performance of the metabolomic PLS-DA classification was compared with that of LVEF and NT-proBNP. The LVEF classification followed the European Society of Cardiology guidelines (34), identifying three classes: reduced LVEF (<40%, high-risk: HiR), mid-range LVEF (40–49%, intermediate risk: IR), and preserved LVEF (≥50%, low risk: LR). Patients with a baseline level of NT-proBNP higher than 400 pg/ml were considered at high-risk of death. Moreover, the ability of the combination of metabolomics and LVEF in predicting poor prognosis was also tested. Metabolomics was used as the first screening method and then was adjusted using ejection fraction: patients predicted as survivors by the PLS-DA metabolomic model but with LVEF of <35% were reclassified as high-risk of death, while patients predicted as deceased but with LVEF of >50% were reclassified as low risk of death. All the above-mentioned analyses were performed using Kaplan–Meier (KM) curves, with the additional calculation of the hazard ratio (HR) and p-value assessed by the Log-Rank test (R library “survminer”). The performances and the independence of metabolomics were evaluated by calculating Cox proportional hazards regression models (41) (R library “Survival”) and each model significance was assessed through a likelihood-ratio test. These analyses were performed in a univariate and multivariate fashion.

The untargeted quantification of 22 metabolites and 114 lipoprotein-related parameters was performed using the Bruker IVDr analysis platform (42). The non-parametric Wilcoxon Rank-Sum test was used to infer differences between the groups of interest. The P-values were adjusted for multiple testing using the false discovery rate (FDR) procedure with Benjamini and Hochberg (43) correction at α = 0.05. Each metabolite/lipoprotein feature was divided into three tertiles and Cox regression models were calculated to estimate the association between metabolites/lipoproteins and prognosis. Additional models were calculated to adjust for additional covariates: sex, age at DCM diagnosis, time from DCM diagnosis and last follow up, NT-proBNP, LVEF, New York Heart Association (NYHA) class at enrollment, systolic blood pressure (SBP), end-diastolic diameter index (EDDi) at enrollment, and left atrial volume index (LAVi) at enrollment.

Robust correlations were calculated among metabolomic variables and clinical data following the 10% winsorized correlation approach (44) using the function “wincor” of the R package “WRC2.” The P-values were adjusted for multiple testing using the FDR procedure.

Baseline characteristics of the cohort are shown in Table 1. The median age at enrollment was 58.5 ± 14.1 years, with a median age at DCM diagnosis of 49 ± 11.9 years. The patients were predominantly men (69.8%). The median time from DCM diagnosis at last evaluation was 15 ± 5.9 years. The patients were mostly pauci-symptomatic at enrollment (85.8% NYHA class I-II). Median systolic blood pressure (SPB) was 120 ± 14.8 mmHg, while echocardiographic parameters showed a considerable LV (left ventricular) and LA (left atrial) enlargement in most patients, with median indexed LVEDDi and LAVi of 31.6 ± 4.9 and 41.9 ± 18 mm/m², respectively. Median LVEF at enrollment was 44.5 ± 8.15% with an NT-proBNP median value of 219.1 ± 237.8 pg/ml. All patients were on ACE-inhibitors or angiotensin receptor blockers inhibitors (ARBs) at enrollment, while 92.4% were receiving beta-blockers (BB). The introduction of Angiotensin Receptor Neprilysin Inhibitors (ARNIs) into clinical practice has become significant only in the last 2 years of the study follow-up, therefore, this agent has not been taken into consideration. Diuretic treatment was administered in 71% of patients at baseline.

At the last evaluation, 80 patients were alive, while 26 had died for HF-related causes (7 patients for sudden cardiac death and 19 patients for refractory heart failure). Demographic and clinical characteristics of the 106 patients with DCM depending on survival status at t2 are shown in Table 1: age, SBP, EDDi, LVEF, NT-proBNP, and diuretic administration at enrollment were significantly different between survivors and deceased patients.

Before any supervised approach, PCA was used as an exploratory analysis to visualize data. No outlier or relevant clustering emerged from the score plot (Supplementary Figure S2). Then, the differences in serum metabolomics ¹H NMR fingerprints of survivors and deceased patients were analyzed using standard PLS-DA (Figure 2A). To ensure that the calculated PLS-DA model was statistically robust, an internal validation using a LOOCV was performed: the two groups show a good clustering, yielding 70.8% accuracy, 76.9% sensitivity, and 68.8% specificity (deceased patients wrongly classified as survivors present, on average, a longer survival time between enrollment and death: 5.7 vs. 4.3 years, see Supplementary Figure S1). Based on the VIP score analysis (Figure 2B) variables that mainly contributed to the PLS-DA model, and, thus, to the risk of stratification, were related to spectral regions of creatine, creatinine, lactate, trimethylamine-N-oxide, and lipoproteins.

Analyzing the metabolomic classification with KM curves (Figure 2C), clear discrimination resulted between patients who died and those who survived: p-value of 0.00002 and HR of 5.71 (95% CI 2.57–12.67). The use of the majority vote algorithm for patient classification produced a gray region in which patients were assigned to a class or another with a small margin (50 ± 10%). To verify the robustness of our model, we repeated the KM analyses after removing the 11 ambiguous patients that were classified with a margin < 10%. Removing borderline patients did not significantly affect the overall model accuracy (Supplementary Figure S3).

The prognostic significance of NT-proBNP natriuretic peptide, examined at enrollment, was analyzed. The results of the KM analyses are shown in Supplementary Figure S4. Compared to metabolomics results, NT-proBNP showed higher specificity (71.2%) in discriminating deceased and survivor patients, although with lower sensitivity (58.3%) and accuracy (68%). The combination of NT-proBNP and metabolomics provided only a slight improvement in the outcome prediction: 75% sensitivity, 68.5% specificity, and 70.1% accuracy.

The predictive performance of the metabolomic model was also compared with that of left ventricular ejection fraction, and the results are reported in Figure 3. Patients with reduced LVEF showed a significantly higher risk of death (HRs of 7.47 and 4.26 and p-values of 0.0001 and 0.0002 concerning LR and IR, respectively), whereas patients with mid-range and preserved LVEF were both associated with a lower risk (HR 1.79, p-value 0.38).

The potential of metabolomics in sub-stratifying LVEF classes was then tested: mid-range and preserved LVEF were considered as a single low-risk class, whereas reduced LVEF was considered as a high-risk class. Although metabolomics showed to effectively sub-stratify low- and high-risk patients in both classes (HRs 3.46 and 6.01, p-values 0.03 and 0.004, respectively), it showed the best performance in the high-risk class (Figure 4).

Finally, the outcome prediction of the combination of metabolomics and LVEF was tested. The resulting combined score relies on LVEF for <35 or >50% values (i.e., when LVEF appears highly prognostic), whilst it is based on metabolomics for intermediate LVEF values, for which LVEF shows the least predictive capacity. This combined approach led to an improved prediction, with 75.5% accuracy, 73.8% sensitivity, 80.8% specificity, HR 8.09, and p-value 0.0000004 as shown by KM analysis (Figure 5).

Possible covariates, such as advanced age at diagnosis (>60 years), sex, the time elapsed since the first diagnosis, NYHA functional class, SBP and echocardiographic parameters (EDDi and LAVi), were compared with the three main prognostic parameters (metabolomics, NT-proBNP, and LVEF) via univariate and multivariate Cox regression analyses (results displayed in Table 2). After univariate analysis, metabolomics, age at diagnosis, advanced NYHA classes (III-IV) at enrollment, NT-proBNP, and reduced LVEF resulted to be statistically associated with the outcome; however, after multivariate Cox regression model, only metabolomics (HR 6.72, p 0.0007), advanced age at diagnosis (HR 3.26, p 0.008), and reduced LVEF (HR 2.6, p 0.001) maintained their correlation with prognosis, regardless of other variables.

Correlations among metabolomics (metabolites and lipoprotein main parameters) and clinical variables are shown in Supplementary Figure S5. Several statistically significant correlations emerged from this analysis. In particular, succinic acid, acetone, and trimethylamine-N-oxide show a strong correlation pattern with NT-proBNP and with several echocardiographic (ECO) parameters (ECO left the atrial end-diastolic area, ECO left atrial end-diastolic diameter, ECO left atrial end-diastolic, and systolic volume). Also, glucose, lactic acid, and citric acid present correlations with the abovementioned echocardiographic parameters. Glycated hemoglobin (Hb1AC) significantly correlates with 12 metabolites. The LVEF anticorrelates with trimethylamine-N-oxide and correlates with creatine, Apo-A1, and Apo-2. Trimethylamine-N-oxide and creatine show the same trend of correlation with systolic blood pressure. Estimated glomerular filtration rate anticorrelates with creatinine and correlates with high-density lipo-protein (HDL) cholesterol, body mass index (BMI) correlates with phenylalanine and tyrosine, and electrocardiogram QT interval shows positive correlations with valine, pyruvic and citric acids.

Several statistically significant correlations also emerged among quantified metabolites and many lipoprotein-related parameters. In particular, amino acids, as well as lipoproteins, are shown to be highly intercorrelated. All results are presented in Supplementary Figure S6.

Univariate analysis of the quantified metabolomic features (Supplementary Table S1), using the Wilcoxon test, unraveled that deceased patients as compared with survivors were characterized at enrollment by lower levels of creatine, apolipoprotein (apo-)A2 HDL, apo-A2, phospholipids HDL-3, apo-A1 HDL-3, apo-A1 HDL-4, apo-A2 HDL-4, phospholipids VLDL, phospholipids HDL-4, apo-A2 HDL-3, free cholesterol VLDL-3, triglycerides HDL-4, cholesterol HDL-3, and by higher levels of trimethylamine-N-oxide, creatinine, lactate, LDL and HDL Cholesterol ratio, triglycerides LDL-3, and triglycerides LDL-2 (p-value <0.05 for all before FDR correction).

Concentrations of each metabolite/lipoprotein were used to build Cox regression models for the evaluation of their net effect on survived and deceased patients (Table 3). Multivariate models adjusted for sex, age at DCM diagnosis, time from DCM diagnosis and last follow up, NT-proBNP, LVEF, NYHA class at enrollment, SBP, EDDi at enrollment, and LAVi at enrollment were also calculated. In the univariate model, higher levels of trimethylamine-N-oxide (3rd tertile), creatinine (3rd tertile), acetic acid (2nd tertile), succinic acid (3rd tertile), triglycerides low-density lipoprotein (LDL) (2nd tertile), and triglycerides LDL-3 (3rd tertile) are associated with a higher risk of poor outcome, whereas higher levels of triglycerides very-low-density lipoprotein (VLDL) (3rd tertile), Apo-A2 HDL (2nd and 3rd tertiles), triglycerides, VLDL-3 (3rd tertile), phospholipids (VLDL-1 (3rd tertile), phospholipids VLDL-3 (3rd tertile), and triglycerides HDL-4 (3rd tertile) are associated with a good prognosis. Among them, only the associations related to trimethylamine-N-oxide (3rd tertile), triglycerides LDL (2nd tertile), phospholipids VLDL-3 (3rd tertile), triglycerides LDL-3 (3rd tertile), and triglycerides HDL-4 (3rd tertile) remain statistically significant in the multivariate model adjusted for clinical covariates. However, in addition to these variables, additional ones (that were not significant in the univariate model) related to various LDL subfractions 1, cholesterol HDL-2, and Apo-B100 Apo-A1 ratio resulted to be significant in the multivariate model (Table 3).