Sample collection

For this study, patients were drawn from two simultaneously recruiting studies in order to capture both stable liver transplant patients and those presenting for a liver biopsy for possible rejection. Both studies were approved by the Emory University Institutional Review Board and in line with all research guidelines. Inclusion criteria were (1) post liver transplant, (2) plasma sample available and (3) either rejection or no rejection by medical record review at the time of the plasma sample collection. We included 43 plasma samples from two studies. Study 1 was the Evaluation of cardiovascular risk markers in pediatric transplant recipients study (PI Vos; IRB#00076255). Patients, ages 10–21 years old, who had undergone a liver transplant at least 12 months prior and who did not have ACR within the last 3 months were included in this study. Study 2 was the Serum markers and magnetic resonance imaging (MRI) in the evaluation of liver disease study (PI Vos, IRB#00002117 and 00094514). All participants were assented, and informed consent was obtained from parents/guardians to participate. This study enrolled any child scheduled for liver biopsy, without fever (in prior 2 weeks) and no chronic renal disease or insufficiency. ACR was determined by pathology and no ACR was determined by medical chart review for any rise in liver enzymes resulting in a clinical rejection diagnosis or liver biopsy within 1 month prior to or after the research visit. Venous blood draws were completed between 08/2014–02/2018 at time of liver biopsy using ethylenediaminetetracetic acid (EDTA) blood tubes for collection of plasma. Demographic questionnaires were also collected at this time. EDTA tubes were inverted several times and immediately put on ice. Plasma was centrifuged for 10–20 min at 1200–1300 rcf, aliquoted and immediately frozen at − 80 °C.

HRM methods

HRM was completed using established methods^25,26. Plasma samples were prepared and analyzed in batches of 20; each batch included duplicate analysis of pooled human plasma (QStd-3) for quality control purposes and reference standardization. Prior to analysis, plasma aliquots were removed from storage at − 80 °C and thawed on ice. Each cryotube is then vortexed briefly to ensure homogeneity, and 50 μL transferred to a clean microfuge tube. Immediately after, plasma was treated with 100 μL of ice-cold LC–MS grade acetonitrile (Sigma Aldrich) containing 2.5 μL of internal standard solution with eight stable isotopic chemicals selected to cover a range of chemical properties. Following addition of acetonitrile, plasma is then equilibrated for 30 min on ice, upon which precipitated proteins are removed by centrifuge (16.1×g at 4 °C for 10 min). The resulting supernatant (100 μL) is removed, added to a low volume autosampler vial and maintained at 4 °C until analysis (< 22 h).

Samples were analyzed in triplicate using 10 μL injections and separate HILIC and C₁₈ chromatography columns with detection by high-resolution mass spectrometry (Q-Exactive HF Orbitrap, Thermo Scientific, San Jose, CA). During HILIC chromatography, the electrospray ionization (ESI) source is operated in positive ion mode while the reverse phase column is flushing with wash solution. Flow rate is maintained at 0.35 mL/min until 1.5 min, increased to 0.4 mL/min at 4 min and held for 1 min. Solvent A is 100% LC–MS grade water, solvent B is 100% LC–MS grade acetonitrile and solvent C is 2% formic acid (v/v) in LC–MS grade water. Initial mobile phase conditions are 22.5% A, 75% B, 2.5% C hold for 1.5 min, with linear gradient to 77.5% A, 20% B, 2.5% C at 4 min, hold for 1 min, resulting in a total analytical run time of 5 min. During the flushing phase, the HILIC column is equilibrated with a wash solution of 77.5% A, 20% B, 2.5% C.

The C₁₈ column is operated parallel to the HILIC column. During operation of the C₁₈ method, the ESI source is operated in negative ion mode while the HILIC column is flushing with wash solution. Flow rate is maintained at 0.4 mL/min until 1.5 min, increased to 0.5 mL/min at 2 min and held for 3 min. Solvent A is 100% LC–MS grade water, solvent B is 100% LC–MS grade acetonitrile and solvent C is 10 mM ammonium acetate in LC–MS grade water. Initial mobile phase conditions are 60% A, 35% B, 5% C hold for 0.5 min, with linear gradient to 0% A, 95% B, 5% C at 1.5 min, hold for 3.5 min, resulting in a total analytical run time of 5 min. During the flushing phase (HILIC analytical separation), the C₁₈ column is equilibrated with a wash solution of 0% A, 95% B, 5% C until 2.5 min, followed by an equilibration solution of 60% A, 35% B, 5% C for 2.5 min.

The high-resolution mass spectrometer was operated in full scan mode at 120,000 resolution and mass-to-charge ratio (m/z) range 85–1275. Probe temperature, capillary temperature, sweep gas and S-Lens RF levels were maintained at 250 °C, 300 °C, 1 arbitrary units (AU), and 45 AU, respectively, for both polarities. Positive tune settings for sheath gas, auxiliary gas, sweep gas and spray voltage setting were 45 AU, 25 AU and 3.5 kV, respectively; negative settings were 45 AU, 5 AU and -4.0 kV. Raw data files were extracted and aligned using apLCMS²⁷ with modifications by xMSanalyzer²⁸. Uniquely detected ions consisted of accurate mass m/z, retention time and ion abundance, referred to as m/z features.

Statistical analysis

Descriptive statistics were used to evaluate demographics using Student’s t test and chi square where appropriate. Metabolites were first filtered based on coefficient of variation (CV) and Pearson correlation between technical replicates. Only features that have a median CV less than 50% and the samples with Pearson correlation greater than 0.7 are used for further analysis. The technical replicates are averaged following the quality assessment and only features with at least 80% signal in either the rejection group or no rejection group were retained. Metabolite data were then log2 transformed and quantile normalized to reduce the effect of technical errors on downstream statistical analysis and interpretation. Hypothesis testing included one-way repeated measures analysis of variance (ANOVA) using LIMMA via the xmsPANDA R package version 1.0.7.4⁶. The p-values were adjusted for multiple comparisons using Benjamin Hochberg false discovery rate (FDR) procedure. A fivefold cross-validation was used in the OPLS-DA analysis and the mean cross-validation accuracy and corresponding standard deviation were reported. Generalized linear regression was used to compare groups and control for days between biopsy and transplant. In order to explore the direct comparison between ACR and no ACR, we used orthogonal partial least squares discriminant analysis (OPLS-DA) and used a Variable Importance in Projection (VIP) > 2 for further annotation. Type 1 (− log₁₀ p vs m/z) and Type 2 (− log₁₀ p vs retention time) Manhattan plots were used to visualize the pattern of differential expression across all features with respect to molecular mass and chemical properties, respectively. We conducted univariate linear regression analyses using PLS component scores as the outcome and phenotypes of interest as the independent variables to determine whether the OPLS-DA analysis was successfully identifying those with rejection and without rejection without influence by other factors.

High resolution metabolomics (95 overlapping significant metabolites), clinical data and patient demographics were integrated using xMWAS package in R²⁹. Clinical and demographic data included age, height, weight, BMI, biological sex, race, rejection status, liver disease diagnosis, liver enzymes, albumin, total bilirubin, hemoglobin, platelet count, days since transplantation, and steroid use. Integrative network analysis was performed using sparse partial least squares regression analysis, a multivariate approach for data integration that included associations with |r| > 0.4 and p-value < 0.05. The multilevel community detection method in xMWAS was used for identifying communities of tightly connected clinical and demographic data and significant metabolites that differentiated rejection.

Metabolite annotation

Metabolic features were annotated using xMSannotator in which the confidence scores for annotation are derived from a multi-stage clustering algorithm³⁰. Further identification of the selected metabolites were confirmed by criteria of Schymanski et al.³¹ either by Level 1 identification, which involves comparing mass spectrum and co-elution relative to authentic standards within a 30-s retention time window, or by Level 2 identification, which involves comparison to METLIN spectral database (http://metlin.scripps.edu/index.php). Lower confidence annotations designated as Level 3–5 identification by Schymanski et al.³¹ were made using HMDB (Human Metabolome Database, http://www.hmdb.ca/)³² and KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.genome.jp/kegg/)³³. Additional manual search was done using METLIN at 5 ppm tolerance³⁴. Only metabolites corresponding to Level 1 identification are reported in this manuscript.

Mummichog v2.0 was used to perform pathway enrichment analysis using m/z features that were significant at p < 0.05³⁵. Mummichog was designed to perform pathway and network analysis for untargeted metabolomics. The software compares the enrichment pattern of the significant metabolite subsets with null distribution on known metabolic reactions and pathways, thereby allowing prioritization of pathways for further evaluation³⁶. Previously published studies have shown that FDR correction results in type 2 statistical error while protecting for type I statistical error²⁵. Pathway enrichment analysis using features significant at raw p-value, provides a 2 step approach which protects against both type I and type II errors³⁶.