Metabolomics biomarkers of hepatocellular carcinoma in a prospective cohort of patients with cirrhosis

Background & Aims The effectiveness of surveillance for hepatocellular carcinoma (HCC) in patients with cirrhosis is limited, due to inadequate risk stratification and suboptimal performance of current screening modalities. Methods We developed a multicenter prospective cohort of patients with cirrhosis undergoing surveillance with MRI and applied global untargeted metabolomics to 612 longitudinal serum samples from 203 patients. Among them, 37 developed HCC during follow-up. Results We identified 150 metabolites with significant abundance changes in samples collected prior to HCC (Cases) compared to samples from patients who did not develop HCC (Controls). Tauro-conjugated bile acids and gamma-glutamyl amino acids were increased, while acyl-cholines and deoxycholate derivatives were decreased. Seven amino acids including serine and alanine had strong associations with HCC risk, while strong protective effects were observed for N-acetylglycine and glycerophosphorylcholine. Machine learning using the 150 metabolites, age, gender, and PNPLA3 and TMS6SF2 single nucleotide polymorphisms, identified 15 variables giving optimal performance. Among them, N-acetylglycine had the highest AUC in discriminating Cases and Controls. When restricting Cases to samples collected within 1 year prior to HCC (Cases-12M), additional metabolites including microbiota-derived metabolites were identified. The combination of the top six variables identified by machine learning (alpha-fetoprotein, 6-bromotryptophan, N-acetylglycine, salicyluric glucuronide, testosterone sulfate and age) had good performance in discriminating Cases-12M from Controls (AUC 0.88, 95% CI 0.83-0.93). Finally, 23 metabolites distinguished Cases with LI-RADS-3 lesions from Controls with LI-RADS-3 lesions, with reduced abundance of acyl-cholines and glycerophosphorylcholine-related lysophospholipids in Cases. Conclusions This study identified N-acetylglycine, amino acids, bile acids and choline-derived metabolites as biomarkers of HCC risk, and microbiota-derived metabolites as contributors to HCC development. Impact and implications: The effectiveness of surveillance for hepatocellular carcinoma (HCC) in patients with cirrhosis is limited. There is an urgent need for improvement in risk stratification and new screening modalities, particularly blood biomarkers. Longitudinal collection of paired blood samples and MRI images from patients with cirrhosis is particularly valuable in assessing how early blood and imaging markers become positive during the period when lesions are observed to obtain a diagnosis of HCC. We generated a multicenter prospective cohort of patients with cirrhosis under surveillance with contrast MRI, applied untargeted metabolomics on 612 serum samples from 203 patients and identified metabolites associated with risk of HCC development. Such biomarkers may significantly improve early-stage HCC detection for patients with cirrhosis undergoing HCC surveillance, a critical step to increasing curative treatment opportunities and reducing mortality.

to upload the related de-identified image files and documentation to a secure dedicated server space.

Global Metabolomics Profiling
Samples were prepared using the automated MicroLab STAR ® system (Hamilton Company).For quality control (QC) purposes, several recovery standards were added prior to the first step in the extraction.
Proteins were precipitated with methanol under vigorous shaking for 2 minutes (Glen Mills GrnoGrinder 2000) followed by centrifugation.The resulting extract was divided into five fractions: two for analysis by two separate reverse phase (RP)/UPLC-MS/MS methods with positive ion mode electrospray ionization (ESI), one for analysis by RP/UPLC-MS/MS with negative ion mode ESI, one for analysis by HILIC/UPLC-MS/MS with negative ion mode ESI, and one sample was reserved for backup.
Samples were briefly placed on a TurboVap® (Zymark) to remove the organic solvent.The sample extracts were stored overnight under nitrogen before preparation for analysis.Several types of controls were analyzed with the experimental samples: a pooled matrix sample generated by taking a small volume of each experimental sample served as a technical replicate throughout the data set; extracted water samples served as process blanks; and a cocktail of QC standards.Instrument variability was determined by calculating the median relative standard deviation (RSD) for the standards that were added to each sample prior to injection into the mass spectrometers.Overall process variability was determined by calculating the median RSD for all endogenous metabolites present in 100% of the pooled matrix samples.
Experimental samples were randomized across the platform run with QC samples spaced evenly among the injections.All methods utilized a Waters ACQUITY ultra-performance liquid chromatography (UPLC) and a Thermo Scientific Q-Exactive high resolution/accurate mass spectrometer interfaced with a heated electrospray ionization (HESI-II) source and Orbitrap mass analyzer operated at 35,000 mass resolution.
The sample extract was dried then reconstituted in solvents compatible to each of the four methods.Each reconstitution solvent contained a series of standards at fixed concentrations to ensure injection and chromatographic consistency.One aliquot was analyzed using acidic positive ion conditions, chromatographically optimized for more hydrophilic compounds.In this method, the extract was gradient eluted from a C18 column (Waters UPLC BEH C18-2.1x100mm, 1.7 µm) using water and methanol, containing 0.05% perfluoropentanoic acid (PFPA) and 0.1% formic acid (FA).Another aliquot was also analyzed using acidic positive ion conditions, however it was chromatographically optimized for more hydrophobic compounds.In this method, the extract was gradient eluted from the same afore mentioned C18 column using methanol, acetonitrile, water, 0.05% PFPA and 0.01% FA and was operated at an overall higher organic content.Another aliquot was analyzed using basic negative ion optimized conditions using a separate dedicated C18 column.The basic extracts were gradient eluted from the column using methanol and water, however with 6.5mM Ammonium Bicarbonate at pH 8.The fourth aliquot was analyzed via negative ionization following elution from a HILIC column (Waters UPLC BEH Amide 2.1x150 mm, 1.7 µm) using a gradient consisting of water and acetonitrile with 10mM Ammonium Formate, pH 10.8.The MS analysis alternated between MS and data-dependent MS n scans using dynamic exclusion.The scan range varied slighted between methods but covered 70-1000 m/z.Raw data was extracted, peak-identified and QC processed using Metabolon's hardware and software.Compounds were identified by comparison to library entries of purified standards or recurrent unknown entities.

Statistical analyses
Table S1: Two-tailed t-test for continuous variables and Fisher test for categorical variables were used to compare demographic and clinical parameters between patients enrolled at the two sites as well as between patients who developed HCC during follow-up (Cases) and those who didn't (Controls).
Tables S2 and S3, Figures 1A and S5, Supplementary Figure and Table: To determine the association between circulating metabolite levels and HCC outcome (at any time point; within 6 months, 12 months or 24 months prior to HCC diagnosis) or HCC treatment, linear mixed-effects modeling using Maximum Likelihood Estimation was performed, using the "lmer" function of the "lme4" R package.Grouping based on binary clinical outcome (Cases versus Controls for HCC development vs no development; Cases-T vs Cases for HCC patients post-versus pre-treatment), time (in months to latest visit) and the time × group interaction term were modeled as fixed effects, whereas patient ID was modeled as a random effect to account for repeated measures from the same patient.Log-transformed metabolite levels were used as the outcome variables.For each metabolite, the model was represented by [y = β0 + β1x1 + β2x2 + β3x1x2 where y is the log-transformed metabolite abundance, β0 is the global intercept, β1 is the coefficient for time (x1), β2 is the coefficient for group (x2) (binary clinical outcome, with Controls=0 and Cases=1; or Cases=0 and Cases-T=1), β3 is the coefficient for time × group interaction (i.e.difference between cases vs controls), µ is the patient-specific random effect, and ϵ is the residual error.The p-value of each fixed effect was calculated using normal approximation and was also corrected for multiple testing by the Benjamini-Hochberg method, giving q-values.The coefficient, p-value and q-value for the fixed effect "group" were used for subsequent volcano plots.

Figures S2 and 1B:
To identify highly correlated clusters of HCC-associated metabolites, Spearman's correlation was performed between the 150 metabolites associated with HCC.

Figure 1C:
Logistic regression was performed to determine the accuracy of each metabolite in predicting HCC development after adjusting for clinical covariates.Using the "glm.fit"function, we obtained odds ratios adjusted for age, gender and diabetes (AOR) and 95% confidence intervals for each unit increase in the level of metabolite (normalized, imputed, log-transformed values).

Figures S3, S6, 2A and 3A:
To determine the predictive ability of a selected combination of HCC-related metabolites, the conditional inference random forest machine learning algorithm was implemented, using the "cforest" function of the "party" package, combined with the "caret" package.Demographic and genetic variables (age, gender, PNPLA3 rs738409 GG genotype, TM6SF2 rs58542926 CT/TT genotype), AFP and the HCC-associated metabolites were used as independent variables, while HCC diagnosis was used as the binary outcome.Using the "train" function of the "caret" package, the optimal "mtry" value giving the maximum AUC was determined by 5-fold cross validation.Conditional importance scores (Mean Decrease in Accuracy) at the optimal "mtry" were generated to obtain final rankings of importance.

Figures 2B-C and 3B-C:
To further select for top variables, recursive feature elimination was implemented using the "rfe" function of the "caret" package and incorporating resampling through 3-fold cross validation.The number of variables giving optimum model performance by the "accuracy" metric was determined.Receiver operating characteristic (ROC) curve analysis was performed to determined the predictive accuracy of a minimal panel, consisting of the 6 predictors with the highest importance by recursive feature elimination.Logistic regression models were fit using the "glm" function in R; the resulting fitted probabilities were used for graphing ROC curves and computing the area under the curve (AUC) using the pROC and ROCR packages.4: To determine the effect of demographic, clinical or genetic variables on HCCassociated metabolomic profiles, redundancy analysis (RDA) was performed using the "capscale" function in the Vegan package for R. Normalized, imputed, log-transformed values of the HCC-associated metabolites were used as the response variables, while age, gender, etiologies, PNPLA3 rs738409 (in increasing number of the risk allele G) and TM6SF2 rs58542926 (in increasing number of the risk allele T), were used as the explanatory variables.To identify the specific metabolites affected by gender, TM6SF2 rs58542926 and PNPLA3 rs738409, linear mixed-effects modeling was again performed.For each of the three clinical variables, the model was represented by [y = β0 + β1x1 + µ + ϵ], where y is the log-transformed metabolite abundance, β0 is the global intercept, β1 is the coefficient for the clinical variable (x1) (gender: female=0, male=1; TM6SF2 rs58542926: CC=0, CT=1, TT=2; PNPLA3 rs738409: CC=0, GG=1), µ is the patient-specific random effect, and ϵ is the residual error.For each model, the coefficient, p-value and Benjamini-Hochberg-adjusted q-values of the clinical variable was generated.5: To determine whether HCC-associated metabolites could distinguish between patients with LI-RADS-3 lesions that developed HCC (Cases-LR3) and those that did not (Controls-LR3), linear mixed-effects modeling was again performed using only Controls and Cases with LI-RADS-3 lesions.

Table S5, Figure
Log-transformed metabolite levels were used as the outcome variable.For each metabolite, the model was represented by [y = β0 + β1x1 + µ + ϵ], where y is the log-transformed metabolite abundance, β0 is the global intercept, β1 is the coefficient for group (x1) (binary clinical outcome, with Controls=0 and Cases=1), µ is the patient-specific random effect, and ϵ is the residual error.The p-value of each fixed effect was calculated using normal approximation and was also corrected for multiple testing by the Benjamini-Hochberg method, giving q-values.The coefficient, p-value and q-value for the fixed effect "group" were used for subsequent volcano plots.Using the 23 significant metabolites from linear mixedeffects modeling, principal component analysis (PCA) was performed using the "cmdscale" function, and Euclidean distances based on log-transformed metabolite levels.Beta dispersion and permutational multivariate analysis of variance (PERMANOVA) tests were performed with the Vegan package.Ellipses were drawn using the standard deviation of point scores.graphs showing for each pathway (Gr1 to Gr9), the percentage of metabolites detected in 100% of the samples or in over 75% of the samples.The majority (58.1%) of the metabolites in the Amino Acids superpathway were detected in all 612 serum samples, with 89.9% of them detected in at least 75% of the samples.Similarly, the majority of the Lipids-related metabolites (50.8%), of the Carbohydrates-related metabolites (55.5%) and of the Nucleotides-related metabolites (55.8%), were detected in all samples, with 82.7%, 81.5% and 90.7% of them detected in at least 75% of the samples, respectively.The 11 energy-related metabolites were also detected in the large majority of samples (at least 99% of the samples).More variation was observed for the detection of the Cofactors and Vitamins-related metabolites (43.5% detected in all samples, 79.5% detected in at least 75% of the samples) and of the Peptides-related metabolites (47.1% detected in all samples, 90.2% detected in at least 75% of the samples).The largest variations were observed as anticipated for the Partially Characterized Molecules, with 24.1% of them detected in all samples and 68.9% in at least 75% of the samples, as well as for Xenobiotics with only 9.3% detected in all samples and 33.0%detected in at least 75% of the samples.A large number of xenobiotics (33.6%) were detected in less than 5% of the samples.Metabolites with a Spearman's correlation coefficient of r>0.9 and p<0.05 are indicated with a white cross.

Fig. S1 .
Fig. S1.Overall distribution of the metabolites measured in the study.(A) Pie chart displaying the

Fig. S2 .
Fig. S2.Correlation matrix of the 151 differential metabolites between cases and controls.

Fig. S3 .
Fig. S3.Importance of individual clinical parameters and metabolites in predicting HCC diagnosis,

Table S2 : Abundance changes in metabolites between controls and cases.
Significance of metabolites abundance between cases and controls was determined by linear mixed model.Super-pathway groups for

Table S3 . Changes in metabolites abundance in cases collected 12 (Cases-12M) or 6 (Cases-6M) months prior to HCC diagnosis compared to controls.
Significance of metabolites abundance changes between controls and cases-12M or cases-6M was determined by linear mixed model.For these

Table S4 . Metabolites significantly affected by PNPLA3 SNP, TM6SF2 SNP or Gender.
Metabolitesignificance was determined with a linear-mixed model.Super-pathway groups are abbreviated as follows:Gr8: Partially Characterized Molecules.Cases (%) and Cases-T (%) reflect the percentage of detection of each metabolite in cases and cases-T samples, respectively.GPA: Glycerophosphatidic acid, CEHC: carboxyethyl-hydroxychroman, logC: group coefficient (log values).The metabolites are ordered from highest to lowest logC.