Discovery of Targets for Immune-Metabolic Antitumor Drugs Identifies Estrogen-Related Receptor Alpha

Abstract

Drugs that kill tumors through multiple mechanisms have the potential for broad clinical benefits. Here, we first developed an in silico multiomics approach (BipotentR) to find cancer cell-specific regulators that simultaneously modulate tumor immunity and another oncogenic pathway and then used it to identify 38 candidate immune-metabolic regulators. We show the tumor activities of these regulators stratify patients with melanoma by their response to anti-PD-1 using machine learning and deep neural approaches, which improve the predictive power of current biomarkers. The topmost identified regulator, ESRRA, is activated in immunotherapy-resistant tumors. Its inhibition killed tumors by suppressing energy metabolism and activating two immune mechanisms: (i) cytokine induction, causing proinflammatory macrophage polarization, and (ii) antigen-presentation stimulation, recruiting CD8+ T cells into tumors. We also demonstrate a wide utility of BipotentR by applying it to angiogenesis and growth suppressor evasion pathways. BipotentR (http://bipotentr.dfci.harvard.edu/) provides a resource for evaluating patient response and discovering drug targets that act simultaneously through multiple mechanisms.

Significance: BipotentR presents resources for evaluating patient response and identifying targets for drugs that can kill tumors through multiple mechanisms concurrently. Inhibition of the topmost candidate target killed tumors by suppressing energy metabolism and effects on two immune mechanisms. This article is highlighted in the In This Issue feature, p. 517.

Figure 1.

Identification of immune–metabolic regulators. A, Overall schematic of regulation and immune modules of BipotentR. The regulator module identifies regulators of an input pathway using ChIP-seq data. The immune module identifies TFCRs that show immunostimulatory or immunosuppressive properties in bulk tumor transcriptomes and are preferentially active in cancer cells (using single-cell tumor transcriptomes). B, Output of BipotentR regulator module. Potential and significance of regulators to bind cis-regulatory elements of genes in four energy metabolism pathways. Each dot indicates a regulator, colored by individual pathways. OXPHOS, oxidative phosphorylation; TCA, tricarboxylic acid cycle. C, The potential of top predicted master regulators to bind energy metabolism genes. Nuclear receptors are displayed in red. D, TFCRs with positive (or negative) associations with proinflammatory signatures are predicted immunostimulators (purple; or immunosuppressors, orange). E, Top TFCRs predicted to be preferentially active in cancer cells (orange; or CD8⁺ T cells, purple) and their differential activity (estimated from single-cell data). F, Output of BipotentR immune module: combined association with proinflammatory signatures (D, estimate from bulk RNA-seq) and differential activity in cancer cells (E, estimate from single-cell data) are displayed for each TFCR. G, Immune–metabolic regulators identified by BipotentR. Energy regulatory potential (estimated by regulator module) and immune-modulatory potential (estimated by immune module) of TFCRs. Highlighted TFCRs are significant and among the top 15% in both modules. Immunostimulators (purple) and immunosuppressors (orange) are colored. H, Validation of BipotentR-identified targets. Effect of knockout (KO) of target identified by BipotentR on T cell–mediated killing of cancer cells.

Figure 2.

ESRRA inhibition activates antitumor immunity in 4T1 mice. P values using Wilcoxon rank-sum test unless stated otherwise. A, Uniform manifold approximation and projection (UMAP) display of scRNA-seq of tumor-infiltrating CD45⁺ cells from ESRRAi- and vehicle-treated mice. NK, natural killer; Treg, regulatory T cell. B, Markers of activated CD8⁺ T cells in genes differentially expressed by ESRRAi in lymphoid cells from treated mice. Significance of up/downregulation of marker sets estimated using permutation tests. C, Fraction of CD8⁺ T cells identified by flow cytometry. D, Tumor volume comparisons between ESRRAi and control. E, Tumor volume comparisons between ESRRAi with and without CD8 antibody. F, Markers of Tregs in genes differentially expressed by ESRRAi in lymphoid cells from treated mice. Significance of up/downregulation of marker sets estimated using permutation tests. G, Fraction of Tregs identified by flow cytometry. H, Densities of macrophage polarization toward M1 (i.e., for each macrophage cell, macrophage polarization = average expression of M1 markers − average expression M2 markers); see also Supplementary Fig. S2O and S2P. I–K, Measurements done after tumors were surgically removed in ESRRAi- or vehicle-treated mice comparing tumor regrowth rate (i.e., 1—relapse rate; I), lung metastasis deposits (J), circulating tumor cells in the blood (K). PBMC, peripheral blood mononuclear cell. *, P < 0.05; **, 0.05 < P < 0.01; ***, P < 0.001; ns, P > 0.05.

Figure 3.

Signaling induced by ESRRAi in vitro.A, Differentially expressed genes (DEG) between ESRRAi and control in the SKBR3 cell line at three time points. Genes were clustered by K-means. B, Pathway enrichment scores corresponding to clusters of DEGs shown in A. GSVA, gene set variation analysis. C, ESRRA knockout potentiates T-cell killing as observed in CRISPR knockout screens in cancer cells cocultured with T cells. The black line represents the relative position of ESRRA knockout among all gene knockouts ranked from most depleted to least depleted. The significance of ESRRA knockout in screens is also displayed. KO, knockout; SKCM, skin cutaneous melanoma. *, P < 0.05; **, 0.05 < P < 0.01; ***, P < 0.001; ns, P > 0.05.

Figure 4.

ESRRA is activated in immunotherapy-resistant tumors and its inhibition does not adversely affect CD8⁺ T cells. P values estimated by the Wilcoxon rank-sum test. A, ESRRA activity in cancer and immune cells from 30 scRNA-seq cohorts. ALL, acute lymphoblastic leukemia; AML, acute myeloid leukemia; BCC, basal cell carcinoma; BRCA, breast cancer; COAD, colon adenocarcinoma; CHOL, cholangiocarcinoma; DC, dendritic cell; EryPro, erythroid progenitor cell; GBM, glioblastoma multiforme; GMP, granulocyte-macrophage progenitor; HNSC, head and neck squamous cell carcinoma; LIHC, liver hepatocellular carcinoma; MB, medulloblastoma; MM, multiple myeloma; NK, natural killer; NSCLC, non–small cell lung cancer; OPC, oligodendrocyte precursor cell; PAAD, pancreatic adenocarcinoma; pDC, plasmacytoid dendritic cell; SKCM, skin cutaneous melanoma; STAD, stomach adenocarcinoma; UVM, uveal melanoma. B, The chromatin accessibility of ESRRA targets in different cell types from scATAC-seq data of a skin cancer cohort. C, Cancer cell ESRRA activity in patient (skin cancer) tumors with pre– and post–anti–PD-1 treatment for responders and nonresponders. D, Enrichment analysis of an ESRRA-regulated gene set in nonresponding CT26 mice after anti–PD-1 treatment.

Figure 5.

Machine learning evaluation of immune–metabolic targets in predicting patient response to anti–PD-1. A, Overview of the estimation of bipotent target activity score in tumors using a linear function (BTAS). B, Kaplan–Meier plots showing PFS and OS differences for patients receiving anti–PD-1 between the low-risk and high-risk groups defined by the median value of BTAS. C, Survival stratification performance of BTAS versus seven other previously published ICB biomarkers for patients with melanoma receiving ICB. Significance of survival was estimated using log-rank P value (−log₁₀). D–K, Kaplan–Meier plots similar to B, with the low-risk and high-risk groups defined by the median value of the signatures. Signatures displayed are TIDE (D; ref. 97), T-cell inflamed (E; ref. 96), immune (F; ref. 90), cytotoxic (G; ref. 91), IFNG (H; ref. 92), melanocytic plasticity (I; ref. 98), immune checkpoint (J), and PD-L1 expression (K; ref. 96). L, Overview of the estimation of BTAS in tumors using a nonlinear function (deepBTAS). M, Kaplan–Meier plots similar to B with the low-risk and high-risk groups defined by the median deepBTAS. N, Performance improvement (OS and PFS) with the addition of deepBTAS to biomarkers (rows). An improvement quantified as (P value of) increase in the likelihood of multivariate model containing a biomarker and deepBTAS over a model containing the biomarker alone. LRT P, P value from likelihood ratio test. O, Kaplan–Meier plots similar to M done separately on tumors with high (left) and low (right) immune signature.

Figure 6.

Evaluation of bipotent regulators of angiogenesis or growth suppressor evasion. A, Effect of knockout of 28 bipotent targets on T cell–mediated killing of cancer cells. B–M, Performance of BTAS and deepBTAS versus other previously published ICB biomarkers. C, Survival stratification performance of BTAS score versus other ICB biomarkers. Significance of survival was estimated using log-rank P value (−log₁₀). B, D, and E, Kaplan–Meier plots show PFS and OS between the four equal quartiles of risk groups divided based on the following signatures: BTAS (B), immune (D; ref. 90), and cytotoxic (E; ref. 91). F–J, Kaplan–Meier plots show PFS and OS between the four equal quartiles of risk groups divided based on the following signatures: IFNG (F; ref. 92), T-cell inflamed (G; ref. 96), melanocytic plasticity (H; ref. 98), PD-L1 (I; ref. 96), and TIDE (J; ref. 97). K and L, Kaplan–Meier plots show PFS and OS between the four equal quartiles of risk groups divided based on the following signatures: immune checkpoint (K) and deepBTAS (L). M, Performance evaluation of ICB biomarkers in terms of AUC.