Circulating Immune Bioenergetic, Metabolic, and Genetic Signatures Predict Melanoma Patients' Response to Anti–PD-1 Immune Checkpoint Blockade

Activation of antitumor T cells with immune checkpoint inhibitors (ICI) has remarkably improved survival in several types of cancer. However, only a small subset of patients respond to ICI therapy. Also, recurrences are observed in a subset of patients who had initially responded, suggesting the development of acquired resistance. Treatment can also be associated with significant toxicity. Biomarkers are needed to select responsive patients and guide the development of approaches to overcome treatment resistance. Several tumor-associated factors have been associated with clinical outcomes, including PD-L1 expression (1), mutational burden (2), mismatch repair deficiency (3), infiltrating lymphocytes (4), and inflammatory gene expression (5). Obtaining tumors for biomarker assessments can be problematic, particularly in longitudinal, on-treatment studies. Tumor biopsies can carry risks, cause treatment delays, and add cost to patient management. Readily accessible peripheral blood has been the focus of several studies. Although the relationship between immune populations in blood and those in a tumor is not established, frequencies of, for example, circulating myeloid-derived suppressor cells, FOXP3+ T regulatory cells, and eosinophils are associated with clinical outcomes in patients treated with ICIs (6, 7). Although these tumor and blood biomarkers can guide treatment decisions, no biomarker is sufficiently predictive of response in individual patients to be clinically useful. Furthermore, most biomarkers are not readily modifiable to improve treatment outcomes.

Complex metabolic processes within mitochondria regulate T-cell activation. Quiescent T cells have low energy demands and can interchangeably use glucose, lipids, and amino acids to generate ATP through oxidative phosphorylation. Activated cytotoxic CD8⁺ T cells require metabolism shift from oxidative phosphorylation to aerobic glycolysis to maintain effector function (8). The tumor microenvironment represses T-cell mitochondrial biogenesis and function. T cells that infiltrate tumors exhibit significant defects in oxidative metabolism (9). It is also postulated that tumor-infiltrating CD8⁺ T cells lose effector function due to cancer cells outcompeting T cells for glucose consumption (10, 11). Whether circulating immune cells in patients with cancer manifest this dysfunction is not established and how circulating immune cell bioenergetics impact ICI response has not been reported. In diseases influenced by oxidative stress, such as type II diabetes, mitochondrial respirometry measurements of blood cells can serve as surrogates of disease activity (12). In this study, we applied blood cell mitochondrial respirometry measurements as surrogates to predict ICI response. We also integrated plasma metabolomics assessment to understand the potential flux of biochemicals in cells from responder and nonresponder patients undergoing ICI therapy. Furthermore, we incorporated transcriptome assessments to identify possible mechanisms that support metabolite consumption as well as the identification of signatures that have the potential to serve as blood markers associated with response.

Patients with stage III or IV melanoma treated with anti–PD-1 antibodies at the Wake Forest Baptist Comprehensive Cancer Center were studied. All patients provided written informed consent for the research approved by the Wake Forest University Health Sciences Institutional Review Board (IRB), according to the ethical standards put forward by the Belmont Report, federal, state, and local regulations, and policies governing human research. The patient diet was not manipulated. Complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD) were evaluated and assigned according to immune RECIST criteria (13). Twenty patients with melanoma who had manifested CR or PR and 20 with PD were studied in bioenergetics and metabolomics analysis.

Samples were taken at fasting, and whole blood was collected into EDTA before the initial ICI infusion on the day of the initial infusion. PBMCs were isolated using Ficoll–Paque density-gradient centrifugation in accordance with the manufacturer's instructions and resuspended in RPMI 1640 medium supplemented with 2% FBS. For cryopreservation, PBMCs were diluted with a 4°C freezing medium containing 10% DMSO and 90% FBS. Freezing medium was added to cell suspension dropwise under continuous mixing. One milliliter aliquots, 5×10⁶ cells/mL in Nalgene cryotubes (Thermo Fisher Scientific), were placed into a Mr. Frosty freezing container (Thermo Fisher Scientific) and kept overnight at −80°C before being transferred into liquid nitrogen for storage. PBMC were thawed for analysis after 6 to 12 months of storage. The cryopreserved PBMC samples were placed in a 37°C water bath and were gently agitated until they had thawed. Thereafter, the cells were resuspended in 10 mL of RPMI 1640 medium with 2% FBS. After centrifugation at 800 × g for 10 minutes, cells were washed in PBS supplemented with 2% FBS by centrifugation for 10 minutes at 400 × g and, finally, resuspended in DMEM with 10% FBS for 1 hour before analysis. Whole blood was centrifuged to separate plasma and frozen at −80°C.

A Seahorse XF96 extracellular flux analyzer (Agilent) was used to measure OCR and ECAR using the Mitochondrial Stress Test Kit. Cells were treated sequentially with oligomycin, FCCP, and rotenone + antimycin A. ECAR is reported as milli-pH (mpH) units per minute, and OCR is reported as pmol/minute. OCR and ECAR are also normalized against cell counts or expressed as a percentage of the baseline oxygen consumption. Reserve capacity, the difference between the maximal respiration and basal OCR; nonmitochondrial respiration, the difference between measurements taken after the addition of oligomycin and antimycin + rotenone; ATP production, the difference between basal OCR and addition of oligomycin; and proton leak, the difference between ATP-linked and basal OCR (after oligomycin).

PBMCs were incubated with tetramethylrhodamine, ethyl ester perchlorate (TMRE) solution (Abcam) for 30 minutes and analyzed in a microplate reader at a 549 nm excitation and 575 nm emission. For electron microscopy, PBMCs were fixed in glutaraldehyde and examined by electron microscopy using a Tecnai F30 (G2 series) 300kV FEG-TEM, Gatan OneView camera. SerialEM software is used for data acquisition and IMOD for data processing. Mitochondria were also visualized by staining with Mitotraker red and nuclei by Hoechst (Thermo Fisher Scientific); images were acquired with Zeiss LSM 880 confocal microscope with Airyscan, Z-stacked, images were analyzed using Aixia software.

Samples were inventoried and immediately stored at −80°C. Each sample received was accessioned into the Metabolon LIMS system and assigned by the LIMS a unique identifier associated with the source identifier only. This identifier was used to track all sample handling, tasks, results, and so on. The LIMS system tracked the samples (and all derived aliquots). All portions of any sample were automatically assigned their unique identifiers by the LIMS when a new task was created; these samples' relationship was also tracked. All samples were maintained at −80°C until processed. Samples were prepared using the automated MicroLab STAR system from Hamilton Company. Several recovery standards were added before the first step in the extraction process for QC purposes. To remove protein, dissociate small molecules bound to protein or trapped in the precipitated protein matrix. To recover chemically diverse metabolites, proteins were precipitated with methanol under vigorous shaking for 2 minutes (Glen Mills GenoGrinder 2000), followed by centrifugation. The resulting extract was divided into five fractions: two for analysis by two separate reverse phase (RP)/ultrahigh performance liquid chromatography-tandem mass spectroscopy (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 placed briefly on a TurboVap (Zymark) to remove the organic solvent. The sample extracts were stored overnight under nitrogen before preparation for analysis.

All methods used a Waters ACQUITY ultra-performance liquid chromatography (UPLC) and a Thermo Scientific Q-Exactive high resolution/accurate mass spectrometer 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.1×100 mm, 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 aforementioned 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 dedicated C18 column. The basic extracts were gradient eluted from the column using methanol and water, however, with 6.5 mmol/L ammonium bicarbonate at pH 8. The fourth aliquot was analyzed via negative ionization following elution from a HILIC column (Waters UPLC BEH Amide 2.1 × 150 mm, 1.7 μm) using a gradient consisting of water acetonitrile with 10 mmol/L ammonium formate, pH 10.8. The MS analysis alternated between MS and data-dependent MSⁿ scans using dynamic exclusion. The scan range varied slightly between methods but covered 70 to 1,000 m/z. Raw data files are archived and extracted as described below.

Raw data were extracted, peak-identified, and QC processed using Metabolon's hardware and software. These systems are built on a web-service platform utilizing Microsoft's .NET technologies, which run on high-performance application servers and fiber-channel storage arrays in clusters to provide active failover and load-balancing. Compounds were identified by comparison with library entries of purified standards or recurrent unknown entities. Metabolon maintains a library based on authenticated standards that contain the retention time/index (RI), mass-to-charge ratio (m/z), and chromatographic data (including MS/MS spectral data) on all molecules present in the library. Furthermore, biochemical identifications are based on three criteria: retention index within a narrow RI window of the proposed identification, accurate mass match to the library ± 10 ppm, and the MS-MS forward and reverse scores between the experimental data and authentic standards. The MS-MS scores are based on comparing the ions present in the experimental spectrum to the ions present in the library spectrum. Although there may be similarities between these molecules based on one of these factors, all three data points can be utilized to distinguish and differentiate biochemicals. More than 3,300 commercially available purified standard compounds have been acquired and registered into LIMS for analysis on all platforms to determine their analytical characteristics. Additional mass spectral entries have been created for structurally unnamed biochemicals, which their recurrent nature has identified (both chromatographic and mass spectral). These compounds can be identified by acquiring a matching purified standard or classical structural analysis.

Peaks were quantified using AUC. For studies spanning multiple days, a data normalization step was performed to correct variation resulting from instrument interday tuning differences. Essentially, each compound was corrected in run-day blocks by registering the medians to equal one (1.00) and normalizing each data point proportionately. For studies that did not require more than one day of analysis, no normalization is necessary, other than data visualization purposes. In certain instances, biochemical data may have been normalized to an additional factor (e.g., cell counts, total protein was determined by Bradford assay, osmolality, etc.) to account for differences in metabolite levels due to differences in the amount of material present in each sample.

Cell suspension from selected patients were collected and tested for viability (targeting >70% viable cells) and processed for cDNA library construction using the 10× Genomics Chromium platform. Barcoded UMI-tagged libraries are paired-end sequenced on an Illumina NextSeq 500 targeting 1,500 cells per sample, a median read depth of 100,000 reads per cell. Raw bcl and fastq data will be demultiplexed, normalized, and postprocessed using Cell Ranger mkfastq pipelines and QC algorithms (10× Genomics). The Cell Ranger mkfastq pipelines (10× Genomics) were used to process QC and align sequencing reads and optimize cell calling. Graph-based or K-means clustering algorithms implemented in Cell Ranger were applied to delineate populations of diverse immune cell lineages whose identities are frequently confirmed by the presence or absence of marker gene combinations informed by antibody-based cell sorting and identification methods or published transcriptional signatures of FACS-sorted immune cells. The transcript read counts were normalized and statistically compared for differential expression between patients using the negative binomial model with independent dispersions and Chi-square analysis based on proportions of cells positive or negative for gene transcripts. Differentially expressed genes (DEG; FDR-corrected, P < 0.05) were analyzed by ingenuity pathway analysis (IPA) and gene set enrichment analysis (GSEA) to identify gene enrichment for biological processes or signaling pathways that underlie response.

Thawed PBMCs were stained with Fixable Viability Stain 575V (BV605) for 15 minutes at room temperature and washed with PBS+1% FCS. Cells were treated with human FcBlock for 10 minutes to limit nonspecific antibody binding, followed by washing. Cells were then stained with the following antibodies, CD127 BB700, CD16 BV650, CD123 BV711, CD69 BV786, CD3 APC-R700, CD19 APC-H7, Glut-14 PE, CD83 PE-CF594, CD56 PE-Cy5, CD64 PE-Cy7, CD4 BUV395, and CD8 BUV496. Cells were acquired using the Fortessa X20 (BD Biosciences), and data were analyzed using DIVA (BD Biosciences) and FCS Express (De Novo Software) programs.

RNA from cells was obtained using the RNeasy method (Qiagen) according to the manufacturer's direction. qRT-PCR was performed on an ABI Prism 7500 Sequence Detection System (Applied Biosystems). Prestandardized primers and TaqMan probes for LDHC mRNA were used (Applied Biosystems). The reverse transcription and PCR were accomplished using a one-step protocol and TaqMan Universal Master Mix (Applied Biosystems). C_t values were determined, and the relative number of copies (RQ) of mRNA was calculated using the ΔΔC_t method.

For metabolomics studies, a t test was ran on the scaled imputed data; all metabolites with a t test P value <0.1 were considered (n = 70). A novel variation of discrete correlate summation (DCS) was performed on the above metabolites; discrete correlate summation clustering given two groups (responders and nonresponders), correlation matrices of 70 metabolites for 20 subjects each, were constructed using Excel data analysis tools. The correlation matrices were transformed into linear probability matrices using Student t-distribution for (n−2) degrees of freedom, returning two tails. One matrix of log correlation ratios (logcr) was calculated (absolute value of the logarithm of the ratio of each comparison's P value) from the above two matrices. The correlation matrix of this logcr matrix was calculated, and each column was summed. The original logcr matrix was sorted in descending order based on these sums, establishing a clear clustering pattern. The upper cluster of this matrix was determined by the summations of all correlation ratios being above Q1 of the total array, left to right. Principal component analysis (PCA) was performed on the DCS metabolites (DCS + all P < 0.05) for 40 subjects (GraphPad Prism).

The Cell Ranger Single Cell Software Suite v.3.1.0 was used to perform sample de-multiplexing, alignment, filtering, and UMI counting (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger). Student t test was performed to compare profiles between responders and nonresponders.

Metabolites and genes were mapped to their respective pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG; ref. 14). Pathways were annotated with KEGG BRITE hierarchies (https://www.kegg.jp/kegg-bin/show_brite?br08901.keg), and the resultant networks were rendered in the Cytoscape network analysis program (15).

The data generated in this study are available upon request from the corresponding author.

Blood from a group of patients with stage III and IV melanoma who responded to ICI treatment with anti–PD-1 antibody (CR and PR) and a group who did not respond (PD) was compared. Patient groups did not vary in gender, age, body mass index (BMI), blood glucose, liver or kidney function, and other clinical parameters (Table 1). The bioenergetics of PBMC collected prior to treatment were assessed using an Agilent Seahorse bioanalyzer, which measures mitochondrial respiration and glycolytic flux under basal conditions and in response to sequential treatment with electron transport chain complex inhibitors and mitochondrial decouplers (Fig. 1A, N = 20/group). These include the oxygen consumption rate (OCR), which reflects oxidative phosphorylation, and extracellular acidification rate (ECAR), reflecting glycolysis. Surrogates of mitochondrial function were normalized to cell number and are calculated based on the bioenergetic profiles (Supplementary Fig. S1A), including basal mitochondrial function, a qualitative assessment of ATP production, maximal respiration, and respiratory capacity. Patients did not show significant differences in basal OCR (Supplementary Fig. S1B).

Compared with nonresponders, responders manifested higher spare respiratory capacity (26.61 ± 4.7 vs. 9.8 ± 2.1, Fig. 1B), maximal respiratory capacity (383.1 ± 28.8s vs. 251.0 ± 20.26, Fig. 1C), and on average, a 2-fold increase in percent respiration over basal (406.4 ± 61.8 vs. 205.5 ± 12.3, Fig. 1D). Flow cytometry data demonstrated no major differences between groups cell viability, indicating these bioenergetics differences among groups are not due to cell fate (Supplementary Fig. S1C). Examination of mitochondrial RNAs by scRNAseq, which can regulate bioenergetics and mitochondrial biogenesis (16), were not significantly changed between responders and nonresponders (Supplementary Figs. S1D–S1N). Furthermore, staining of mitochondria with mitotracker red in a subset of patients did not show significant differences between patient groups, thus suggesting that changes in mitochondrial bioenergetics may not be due to differences in mitochondrial mass (Supplementary Figs. S1O and S1P). Although we did not observe changes in mitotracker red levels, we observed changes in mitochondrial morphology when examining cells under confocal microscopy (Supplementary Fig. S1O). Mitochondrial dynamics and membrane potential drive metabolic processes, and bioenergetic differences can be associated with mitochondrial polarization and morphology changes. Cristae morphology, in particular, is reflective of the metabolic reprogramming of T cells (17). We examined mitochondrial membrane polarization by treating PBMCs with TMRE. Our data in a subset of patients shows that responder patients have elevated levels of TMRE relative to nonresponder patients (56.6 ± 8.2 vs. 22.6 ± 3.6, Fig. 1E). Depolarized mitochondria are thought to be dysfunctional and fail to uptake TMRE; therefore, these data suggest that the observed OCR changes may be due to defects in mitochondria in nonresponder patients. We also examined patient PBMCs by electron microscopy (EM) and three-dimensional electron tomography (3D-TEM). Mitochondria of responders cells in two-dimensional EM are more spherical; 3D-TEM shows that these mitochondria maintain a reticular shape with “stacked,” well-differentiated cristae; on the other hand, nonresponders also showed a reticular shape, but 3D-TEM shows disorganized expanded cristae (Supplementary Fig. S1Q). Measurement of cristae diameter by 3D-TEM was also observed to be tighter by 3D-EM in responders compared with nonresponders (16.05 ± 0.32 vs. 26.80 ± 0.72, Supplementary Fig. S1R). Circular “fissed” mitochondria are associated with effector cells and increased aerobic glycolysis. Consistent with this, basal ECAR of responders PBMCs were characterized by higher levels of glycolytic flux (8.481 ± 3.4 vs. 5.847 ± 3.3, Fig. 1F), which was significantly correlated with increased overall survival of patients (20.56 ± 8.0 vs. 14.43 ± 9.4, Supplementary Fig. S1S). These bioenergetics parameter changes were associated with a 50% higher ATP production based on oxygen consumption (6.0 ± 3.4 vs. 3.3 ± 2.3, Fig. 1G). These results suggest that response to ICIs may depend on mitochondrial fitness and glycolytic function. Patients that respond to ICI are capable of increasing cellular bioenergetics. Nonresponders may have bioenergetic deficiencies due to damage or anomalies in mitochondrial dynamics.

The cellular bioenergetic results observed by our respirometry analysis could be influenced by the consumption of circulating metabolites in the blood. To determine whether blood metabolites are different between responders and nonresponders, pretreatment plasma was subjected to LC/MS and GC-MS analysis. A total of 1,095 metabolites were identified, with unidentified metabolites in 1.3% and 1.5% in responders and nonresponders, respectively. A t test on the scaled imputed data was performed, and metabolites with a t test P value ≤ 0.1 were considered different (Supplementary Fig. S2). A variation of discrete correlate summation analysis was performed on these metabolites (18), and correlation matrices were transformed into linear probability matrices using Student t-distribution for (n−2) degrees of freedom, returning two tails (Supplementary Fig. S2A). PCA demonstrated that responders cluster together regardless of gender, whereas non-responders segregate from one another (Fig. 2A), suggesting that patients that respond may have similar metabolic profiles. The scaled intensities of each metabolite were compared. We observed, a significant proportion of the metabolites identified derived from lipid, amino acid, carbohydrate, and nucleotide metabolism, while a minor proportion derived from xenobiotics, vitamins/cofactors, drugs, and peptides (Fig. 2B). We identified sub-pathways that differed between metabolites in these pathways that were significantly higher or lower, depending on response, as summarized in Fig. 2C. One of the most significantly upregulated metabolites we observed in plasma was metformin (Fig. 2C; Supplementary Fig. S2B). Further examination showed that 4 responder patients had elevated levels of this drug compared to the rest (Supplementary Fig. S2B). Metformin is implicated in the metabolic reprogramming of immune cells and is implicated in the clinical response to immune checkpoint therapy, including anti–PD-1 therapy (19–21). To determine that this drug metabolite was not driving the observed clustering, we re-analyzed our data without considering this variable. Our PCA data shows that responder and non-responder patients cluster differently even when metformin is not considered a factor (Supplementary Fig. S2C). Therefore, these data suggest that circulating plasma metabolites may be associated with response and may influence cellular energetics of circulating PBMCs.

Lipid metabolites made up 34% of the differentially expressed metabolites identified (Fig. 2C). Long-chain fatty acids scaled metabolic intensities were elevated in responders (P < 0.05, Fig. 3A), with general relative increases including n3 long-chain fatty acids (Fig. 3A and B). Although higher levels of, for example, adrenate (22:2n6), were observed in responders, this metabolite tended to be elevated in both responders and nonresponders when considered within the context of the other metabolites in this pathway (Fig. 3B). Similarly, docosatrienoate (n6 DPA; 22:3n6) was higher in responders than nonresponders but tended to be lower in both groups (Fig. 3B). Nonresponders were characterized by higher levels of myelin-related metabolites (P < 0.05, Fig. 3C), including sphingomyelins (Fig. 3C and D). Significantly higher levels of sphingomyelin-d18:2/24:2 were observed in nonresponders compared with responders (1.3 ± 0.6 vs. 0.96 ± 0.3, P < 0.05, Fig. 3G). Another set of lipids differentially expressed were plasmalogens (P < 0.05, Fig. 3E). There were increases in specific metabolites in this pathway in responders when compared with nonresponders patients (Fig. 3F), including 1-(1-enyl-palmitoyl)-2-arachidonoyl-GPC (P-16:0/20:4; 1.1 ± 0.3 vs. 0.90 ± 0.2, P < 0.01) and 1-(1-enyl-stearoyl)-2-oleoyl-GPE (P-18:0/18:1; 1.2 ± 0.3 vs. 0.97 ± 0.3, P < 0.03; Fig. 3H and I).

Branched-chain amino acids (BCAA) can support mitochondrial function and metabolite consumption differences, explaining the bioenergetics differences observed in patients. BCAA levels were higher in responders than nonresponders (Fig. 4A and B). In general, BCAAs were depleted in nonresponders with no overlap in log-transformed intensities (Fig. 4B). On average there was a 31% increase in isovalerate (i5:0; P < 0.05) and a 35% increase in 2-methyl-butyryl carnitine (P < 0.05) in responders (Fig. 4C and D). BCAAs, which are catabolized in the mitochondria, can regulate inflammatory responses (22). Tryptophan metabolites were also higher in responders (Fig. 4E), although the relative abundance of specific metabolites did vary in responders than nonresponders (Fig. 4F). Tryptophan metabolites, for example, indoleamine 2,3-dioxygenase (IDO), are major regulators of innate and adaptive immune responses (23). Although kynurenine, a marker of IDO activity, was not differentially expressed, xanthurenate, the breakdown product of kynurenine, was higher in responders (Fig. 4G). Other amino acids were significantly regulated between responders and nonresponders. Responders had 38% lower pyroglutamine (Fig. 4H), 45% higher forminoglutamate (Fig. 4I), and 2-fold higher lanthionine (Fig. 4J).

Our PBMC cellular bioenergetics data suggest that ECAR, a surrogate measurement of glycolysis, is elevated in responders versus nonresponder patients (Fig. 1F). Consistent with this finding, levels of glycolytic and pentose phosphate pathway metabolites differed in responders and nonresponders (P < 0.05, Fig. 5A). Nonresponders were characterized by increased pentose phosphate pathway metabolites (Fig. 5A and B). Intermediaries of the pentose phosphate pathway, such as sedoheptulose and ribulonate/xylulonate/lyxonate, were increased in nonresponders by 50% (P < 0.03) and 44% (P < 0.05), respectively (Fig. 5C and D), suggesting differences in sugar transport and metabolism. Although glucose levels were similar (Fig. 5A and B), responders were characterized by increases in glycolytic metabolites, such as 3-phosphoglycerate (Fig. 5A and B). Although pyruvate levels were lower in responders (P < 0.05, Fig. 5E), levels of lactate were increased in responder patients (P < 0.05, Fig. 5F). The ratio of lactate to pyruvate, which can be an indicator of glycolytic activity, was significantly elevated in responder patients (P < 0.05, Fig. 5G), suggesting a predominant glycolytic pattern. The lactate-to-pyruvate ratio was significantly correlated with the overall survival of patients (20.56 ± 8.0 vs. 14.43 ± 9.4, Supplementary Fig. S2D). Because metabolic reprogramming for T-cell activation requires a shift from oxidative phosphorylation to glycolysis, these results indicate that responders may be more primed to execute glycolysis to generate antitumor effectors in response to ICIs.

Single-cell RNA sequencing (scRNAseq) of PBMCs from responder and nonresponder patient pretreatment was performed to examine the molecular regulators of the bioenergetic and metabolomic differences observed. Markers of specific immune cell types, including T cells, natural killer cells, monocytes, dendritic cells (DC), and B cells, were used to annotate immune cell populations among the PBMCs. Significant differences in these immune cell populations between responders and nonresponders were not apparent (Fig. 6A; Supplementary Figs. S3A–S3G). Because the bonafide target of anti–PD-1 are T cells, we examined if there were differences in gene expression within this population associated with response. As portrayed in the volcano plot, responders and nonresponders manifested differential PBMC gene expression changes (Fig. 6B). Consistent with our metabolomic analysis, the most significantly upregulated gene in responders' T-cell population was SLC2A14, which encodes the gene for glucose transporter 14 (Glut-14). T-cell SLC2A14 was upregulated in responders 4-fold compared with nonresponders (P = 2.48e−16, Fig. 6B and C). We then performed flow cytometry in a secondary validation cohort (Table 2), combining a larger subset of patients from those included in the scRNAseq analysis to validate cell populations as identified by scRNAseq (Supplementary Figs. S3H–S3O) and confirm expression of Glut-14 levels. There were no major differences in cell proportions as examined by flow cytometry with the exception of a higher number of CD8⁺ T cells in nonresponders (Supplementary Fig. S3J) and higher CD4⁺ in responder patients (Supplementary Fig. S3K). Our flow cytometry data confirmed the scRNAseq results showing that responders significantly increased Glut-14 in CD3⁺, CD8⁺, and CD4⁺ T cells compared with nonresponders (Fig. 6D–F). In contrast, no significant differences were found in Glut-14 expression in B cells, NK cells, or DC populations (Supplementary Figs. S4A–S4D). Because the cellular bioenergetic studies showed that responders have increases in ECAR (Fig. 1F), a measurement of glycolysis and the scRNA-seq studies showed increased glucose transporter levels, additional genes that regulate glycolytic activity were interrogated. Not surprisingly, LDHC, a gene that encodes lactate dehydrogenase C involved in the conversion of L-lactate to pyruvate, was significantly upregulated in responders' T cells (P < 0.001, Fig. 6G and H). Further examination of PBMC by qRT-PCR in another subset of patients indicated that LDHC was upregulated in responders by almost 50-fold compared with nonresponders (P < 0.001, Fig. 6H). Together with our metabolic analysis, these results suggest that ICI response is associated with a glycolytic signature.

Metabolic abnormalities in the tumor microenvironment limit the ability to mount an antitumor response (10). We found that functional and molecular metabolic biomarkers associated with ICI response can be detected in blood pretreatment. PBMC of patients that respond to treatment have increased ECAR, a measure of glycolytic flux, when compared with patients that do not respond, suggesting that responders have a preference to drive cellular bioenergetics through glycolysis. Patients who respond to therapy also have a common metabolomic signature in their plasma pretreatment, including increased lactate levels to pyruvate and specific lipid and amino acid metabolites. Differences in molecular regulators were also found on single cell transcriptome analysis of circulating PBMCs. The most significant gene upregulated in the T-cell population from responders compared with nonresponders was SLC2A14, which encodes a gene for the glucose transporter, Glut-14. Downstream effectors of the glycolytic pathway, such as LDHC, were also upregulated in responders. Circulating lymphocyte cell mitochondrial morphology also varied. Responders and nonresponders had similar frequencies of circulating immune cell populations on transcriptome and flow cytometric analyses. These suggest that the metabolic differences observed are not due to changes in the proportion of circulating immune cell types but rather to changes in their metabolism.

Circulating immune-cell bioenergetics can differentiate between healthy individuals and those with autoimmune diseases (12, 24, 25). T cells from patients with systemic lupus erythematosus (SLE) are characterized by increases in both oxidative phosphorylation and aerobic glycolysis (26). In the context of SLE, this may represent the chronic activation of immune cells associated with pathogenesis. In the context of ICI, these metabolic changes may be indicative of immune cell activation and antitumor activity. Our data support this possibility, as we observed increased ECAR from circulating PBMCs and SLC2A14/Glut-14, LDHC expression by scRNAseq in the T-cell populations in circulating cells. Although not all plasma metabolites are derived from circulating PBMCs, we observed increased lactate to pyruvate ratio in plasma, suggesting the predominant role of glycolysis in supporting metabolic reprogramming needed to support ICI response. Although baseline OCR was not significantly different, the blockade of ATPase with oligomycin to measure maximal respiration and respiratory capacity was significantly higher in responders. These suggest that the responder immune cells may undertake more stress or engage in metabolic reprogramming more efficiently. Lactate may be a primary carbohydrate fuel in responders. The reverse reaction of lactate-to-pyruvate produces NADH, which can then be transported by the malate-aspartate or the glycerol phosphate shuttles to the mitochondria to support electron transport and, ultimately, ATP production. The lactate-to-pyruvate ratio also functions as a redox buffering system (27); thus, responders patient cells may have the ability to support mitochondrial function under stress by reprogramming the metabolism of fatty acids and amino acids. We observed an abundance of lipid metabolites in responders, including long-chain fatty acids. Serum long-chain fatty acids were associated with ICI response in patients with urologic cancers (28). In this study, enhanced peroxisome signaling in T cells was implicated in reprogramming fatty acid metabolism to overcome metabolic restrictions in the tumor microenvironment. Activation of fatty acid oxidation is implicated in the metabolic reprogramming needed for the development of memory T cells, which are critical to the generation of durable responses to ICIs (29–31). Responders may be primed to support memory T-cell differentiation as they have more abundance of lipid metabolites to support processes, such as β-oxidation. Further studies using samples obtained from patients after cycles of therapy could answer this question.

Lower levels of other lipid metabolites were observed in responders, including several myelin metabolites. This did not include sphingosine-1-phosphate, which in tumors has been associated with ICI resistance (32). Several plasmalogens were increased in responders. These membrane components, which also play roles in cellular signaling and inflammation, are associated with ICI response have not been reported (33). BCAA and glutamate metabolites, which coordinate with fatty acids to modulate energy metabolism and inflammation, were increased in responders (34). Levels of the tryptophan metabolite, kynurenine, increased in patients with melanoma and kidney cancers responding to ICI (35, 36). Although kynurenine was not differentially expressed in our studies, xanthurenate, the breakdown product of kynurenine, was higher in responders.

To maintain effector function, T cells require maintenance of mitochondrial turnover and biogenesis (8, 37). Our previous work in T cells suggested that increases in mitochondrial function are necessary to preserve lymphocyte viability under stress (38, 39). Mitochondrial dynamics are thought to control immune cell metabolism and differentiation (14). Morphologic changes in SLE T-cell mitochondria, as assessed by confocal microscopy and two-dimensional electron microscopy (EM), include increased mitochondrial fission, which functions to mitigate stress, and decreased mitochondrial fusion acts to create new mitochondria (40). Fissed mitochondria with loose cristae are associated with glycolysis, as this morphology would make the electron transport chain less efficient (14). Fused mitochondria with tight cristae are thought to drive oxidative phosphorylation. Responders in our study had lower mitochondrial cristae diameter but were observed to be more spherical as assessed with 3-dimensional EM; the latter is thought to favor glycolysis. Furthermore, the mitochondria in responders showed “stacked,” well-differentiated cristae, suggesting that cristae remodeling more than mitochondrial shape drives a particular metabolic pathway and, in the context of our study, another factor that can predict patient response.

Our study has some limitations. Our scRNAseq data focus on T-cell genes since anti–PD-1 therapy mainly targets T lymphocytes; however, the bioenergetic status and differentiation of other immune cell types could contribute to the retrospective PD-1 response. Still, enrichment analysis of both plasma metabolomics and scRNAseq data suggest regulation of interferon-regulated pathways as well as signatures regulating antigen presentation (Supplementary Fig. S5). The results of this study identify several cellular, soluble, molecular, and metabolomic markers in blood pretreatment that distinguish individuals who will respond to ICI treatment versus those who will not. An ideal biomarker would have both predictive utility and the ability to be targeted therapeutically. Preclinical and clinical studies have demonstrated how several modulators of metabolism can regulate immune responses. Thus, we have provided a comprehensive approach to support personalized treatment strategies that can help guide treatment decisions and develop approaches to address ICI treatment resistance.

E.R. Stirling reports grants from National Institute of General Medical Sciences (NIGMS) during the conduct of the study. Q. Song reports grants from NCI Cancer Center Support Grant during the conduct of the study. W. Zhang reports personal fees from Finnegan, Henderson, Farabow, Garrett & Dunner, and LLP outside the submitted work. D.R. Soto-Pantoja reports grants from The V Foundation for Cancer Research and NIH-NCI during the conduct of the study. No disclosures were reported by the other authors.

P.L. Triozzi: Resources, supervision, investigation, methodology, writing–original draft. E.R. Stirling: Formal analysis, investigation, writing–review and editing. Q. Song: Data curation, formal analysis, validation. B. Westwood: Data curation, software, formal analysis. M. Kooshki: Resources, investigation, methodology. M.E. Forbes: Resources, data curation, formal analysis, methodology. B.C. Holbrook: Software, formal analysis, methodology. K.L. Cook: Formal analysis, validation, investigation, writing–review and editing. M.A. Alexander-Miller: Formal analysis, investigation, writing–review and editing. L.D. Miller: Resources, formal analysis, funding acquisition, validation, visualization, writing–review and editing. W. Zhang: Resources, supervision, funding acquisition, investigation, project administration, writing–review and editing. D.R. Soto-Pantoja: Conceptualization, data curation, supervision, funding acquisition, investigation, writing–original draft, project administration, writing–review and editing.

This work was supported by the V Foundation V Scholar in Cancer Research Award #V2019-018 to D.R. Soto-Pantoja and Wake Forest Baptist Comprehensive Cancer Center Pilot Funds (NIH/NCI P30CA012197) and an NIGMS NRSA T32 Fellowship T32GM127261 (to E.R. Stirling). D.R. Soto-Pantoja was also supported by the American Cancer Society Research Scholar Grant (133727-RSG-19-150-01-LIB) and an NIH-NCI R21 (R21CA249349). The authors wish to acknowledge the support of the Wake Forest Baptist Comprehensive Cancer Center Bioinformatics Shared Resource, Cancer Genomics Shared Resource, Flow Cytometry Shared Resource supported by the NCI's Cancer Center Support Grant award number P30CA012197. We acknowledge the technical contribution of Kenneth Grant of the Wake Forest Cellular Imaging Shared Resource, Dr. Heather Brown-Hardin of the Wake Forest Biology Department Microscopic Imaging Core, and Gary Morgan and Dr. Eileen O' Toole at the University of Colorado, Boulder EM Services Core Facility in the MCDB Department. W. Zhang was supported by the Hanes and Wills Professorship in Cancer and a Fellowship from the National Foundation for Cancer Research. The authors thank patients and families who agreed to participate in the study. The content is solely the authors' responsibility and does not necessarily represent the official views of the NCI.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.