Introduction

Immune checkpoint inhibitors (ICI) antagonizing the programmed cell death 1 (PD-1) signaling axis by disrupting the interaction between PD-1 on exhausted T cells and programmed cell death ligand 1 on tumor and stromal cells have vastly improved outcomes for patients with high-risk and late-stage melanoma. However, among responders to PD-1 blockade therapy, approximately 25% will relapse within 2 years of initiating treatment.1 Little can be done at present for these patients. Analysis of relapsing tumors has revealed that 30% are characterized by intrinsic—that is, ligand-independent—interferon-gamma (IFNγ) signaling.2 Intrinsic IFNγ signaling in cancer cells has been shown to follow chronic IFNγ stimulation and to entail epigenetic remodeling where signal transducer and activator of transcription 1 (STAT1) and interferon regulatory factor 3 (IRF3) permanently occupy the regulatory sites of a subset of interferon-stimulated genes (ISG), including notably those promoting T cell dysfunction.3 Indeed, the ISG selectively upregulated in chronically stimulated cancer cells (ISG resistance signature or ISG.RS) predicts poor survival in patients and poor response to PD-1 blockade.4

Our previous work revealed a contribution from Poly ADP ribosyl polymerase 14 (PARP14) in mediating chronic IFNγ stimulation-driven PD-1 blockade resistance in preclinical models.5 One model we developed entailed implantation of the mouse rapidly accelerated fibrosarcoma homolog B (BRAF) and phosphatase and tensin homolog (PTEN) mutant melanoma cell line YUMM2.16 into the hypodermis of syngeneic C57/BL6 mice, and then treating these mice with a short course of anti-PD-1 (α-PD-1) antibodies when tumors became palpable. In tumors that resumed growth following α-PD-1 treatment, we observed increased expression of multiple ISG, including ISG.RS genes, relative to control (untreated) tumors and to tumors still responding to PD-1 blockade therapy.5 Artificial Intelligence (AI)-assisted immune profiling, based on deconvolving gene expression, indicated that relapsing tumors were also relatively more T cell infiltrated, including by effector memory (TEM) and central memory (TCM) T cells. However, these T cells were also more dysfunctional.5 The regrowth of such tumors was significantly delayed—but not abolished—and survival lengthened by adding a highly potent and selective PARP14 inhibitor RBN012759 (PARP14 catalytic inhibitor (PARP14i), Ribon Therapeutics)7 to the treatment regimen.5

This previous study also highlighted that the PARP14 gene is a STAT1 target and that PARP14 itself is a STAT1 interacting protein.5 Moreover, in macrophages, PARP14 interaction with STAT1 inhibits the expression of proinflammatory genes8 but regulates IRF3 recruitment of RNA polII to promoters of IRF3 target genes9 and is thus positioned to influence the expression of ISG. However, the details are far from clear. Certainly, our previous analysis of gene expression indicated that PARP14 inhibitor treatment of melanoma cells chronically stimulated with IFNγ increased IFN-driven, tumor necrosis factor (TNF)-α-driven, and interleukin (IL)-6-driven inflammatory signaling.5

Intrigued by the potential of PARP14 to influence the expression of ISG in tumor cells but also the failure of PARP14 inhibitor treatment to permanently flatten tumor rebound after α-PD-1 treatment, we undertook a detailed analysis of tumor composition and gene expression changes attending PARP14 inhibitor treatment to identify other targets for intervention which we now present below.

Materials and methodsReagents

Recombinant human IFNγ (PHC4031, Gibco) and recombinant mouse IFNγ (PMC4031, Gibco). PARP14 inhibitor RBN0127597 and PARP14 proteolysis targeting chimera (PROTAC) RBN012811 (Ribon Therapeutics).10

Cell lines

Human melanoma cell lines A375 and 501-Mel (provided by Claudia Wellbrock, The University of Manchester), MC38 cells (provided by Santiago Zelenay, The University of Manchester) and YUMM2.1 cells (provided by Richard Marais, The University of Manchester) were maintained in Roswell Park Memorial Institute (RPMI)-1640 (Sigma-Aldrich) supplemented with 10% v/v fetal bovine serum (Life Technologies) and 1% w/v penicillin-streptomycin (Sigma-Aldrich), at 37°C in a 5% v/v CO2 humidified incubator. Cell lines authenticated by short tandem repat (STR) profiling; cultures routinely tested for mycoplasma. Cells were treated for 2 weeks with 50 IU/mL IFNγ, refreshed every 2–3 days, and subsequently treated with 500 nM Dimethyl sulfoxide (DMSO) (Sigma-Aldrich); 500 nM PARP14i; or 100 nM PARP14 PROTAC for 48 hours.

Mouse tumor implants

Mice were housed in the Biological Services Facility of The University of Manchester on a 12/12 hour light/dark cycle and given ad libitum food (Bekay, B&K Universal, Hull, UK) and water.

8–12 week-old C57BL/6 female mice (ENVIGO) were allowed 1 week to acclimatize. YUMM2.1 cells (5×106 cells) in 100 µL serum-free RPMI-1640 were subcutaneously injected into the left flank. Tumor volumes (height × width × length caliper measurements) and mouse weights were monitored every 2–3 days. When tumors reached a volume of 80 mm3, mice were administered two doses of 300 µg α-PD-1 antibody (BioXCell) in 100 µL InVivoPure pH 7.0 Dilution Buffer (BioXCell) via intraperitoneal (i.p.) injection administered at 3-day intervals. Mice were culled if tumors remained under 100 mm3 (total responders) or conversely reached 800 mm3 (non-responders). Otherwise, mice with tumors recovering to 140 mm3 were randomized to a second round of treatments (retreatment): 2 doses of 300 µg of rat isotype control antibody immunoglobulin G2b (IgG2b) (BioXCell) (IgG group) or α-PD-1 antibody (BioXCell) (α-PD-1 group) in 100 µL InVivoPure pH 7.0 Dilution Buffer (BioXCell) via i.p. injection at 3 day intervals, 14 doses of 500 mg/Kg of PARP14I by oral gavage two times a day (PARP14i group), or a combination of α-PD-1 i.p. and PARP14i gavage in the same doses (Combination group). PARP14i was dissolved in 0.5% w/v methylcellulose (Sigma-Aldrich) + 0.2% v/v Tween 80 (Sigma-Aldrich). Dose volume—10 mL/kg.

Tumor immune infiltrate analysis by flow cytometry

Tumors were dissected on reaching 250–400 mm3 following retreatment and incubated for 45 min with 100 µg/mL Liberase (Sigma-Aldrich) in serum-free media at 37°C, then pushed through a 100 µM nylon cell strainer. The cell suspension was stained for 20 min with Zombie UV Fixable Viability Kit (BioLegend) in phosphate-buffered saline. Subsequently, Fc receptors were blocked, and cells were stained with a surface stain antibody mix. Cells were fixed and permeabilized using the Foxp3/Transcription Factor Staining Buffer Set (eBioscience) following manufacturer instructions. Samples were measured on a BD Fortessa flow cytometer (BD Biosciences), and data was collected using BD FACSDiva software. For all antibodies, a non-stained cell sample and appropriate fluorescence minus one control were analyzed as well. Data were analyzed using FlowJo v8.7.

Antibodies: CD16/32 (clone 93), CD45 (Cat. 103132), CD3ε (Cat. 100306), CD4 (Cat. 100552), CD8α (Cat. 100742), PD-1 (Cat. 135231), LAG-3 (Cat. 125224), TIM-3 (Cat. 119721), Ki-67 (Cat. 652420), Granzyme B (Cat. 372214), CD69 (Cat. 104531), and CTLA-4 (Cat. 106338) from BioLegend. Tox (Cat. 12-6502-82) was purchased from eBioscience. The Gating strategy is described in online supplemental figure 1.

Single-cell RNAseq sequencingTumor processing

Tumors were collected on day 7 of retreatment into serum-free RPMI on ice, then minced into 1 mm pieces and incubated at 37°C up to 60 min in Liberase TL (Millipore Sigma) with rocking according to manufacturer instructions. Digestion was stopped using fetal bovine serum, and tissue was further dissociated using a wide-bore pipette and filtered through a cell strainer to achieve a single-cell solution. Immune population enrichment was done by CD45 mouse Microbeads (Miltenyi) according to manufacturer instructions. Following sorting, both populations were counted, and ∼5% of CD45 negative cells were spiked back into the CD45 positive fraction. Cells were fixed for 18 hours, then quenched and stored at −80°C in Enhancer according to 10x Genomics Tissue Fixation and Dissociation for Chromium Fixed RNA Profiling protocol (CG000553 Rev A) until all samples were collected.

Single-cell isolation and library construction

Gene expression libraries were prepared from formaldehyde-fixed samples using the Chromium X and Chromium Fixed RNA kit, Mouse Transcriptome, 4 rxns × 16 BC (10x Genomics) according to the manufacturer’s protocol (CG000527 Rev C). Illumina-compatible sequencing libraries were constructed by adding P5, P7, i5, and i7 sample indexes and Illumina Small Read 2 sequences to the 10x barcoded, ligated probe products via Sample Index PCR followed by a SPRIselect size-selection.

Sequencing

The resulting sequencing library comprised standard Illumina paired-end constructs flanked with P5 and P7 sequences. The 16 bp 10x Barcode and 12 bp Unique Molecular Identifier (UMI) were encoded in Read 1, while Read 2 sequenced the ligated probe insert, constant sequence, and the 8 bp Probe Barcode. Sample indexes were incorporated as the i5 and i7 index reads. Paired-end sequencing (28:90) was performed on the Illumina NovaSeq 6000 platform.

Data analysisRaw data processing and cell filtering

Raw sequencing data were processed using the 10x Genomics Cell Ranger pipeline (v7.0.0). Base call (BCL) files generated by the sequencer were demultiplexed and converted to FASTQ files using “cellranger mkfastq”. The FASTQ files were then mapped against the pre-built Mouse reference package from 10x Genomics (mm10-2020-A) using “cellranger multi” to demultiplex and produce the gene-cell barcode matr-ix for individual samples. The single-cell data were processed in the R environment (v4.1) following the workflow documented in Orchestrating Single-Cell Analysis with Bioconductor.11 Briefly, for each sample, the Hierarchical Data Format (HDF)5 file generated by Cell Ranger was imported into R to create a SingleCellExperiment object. A combination of median absolute deviation, as implemented by the “isOutlier” function in the scuttle R package (v1.4.0), and exact thresholds were used to identify and subsequently remove low-quality cells before data integration.

Data integration, visualization, and cell clustering

The log-normalized expression values of the combined data were re-computed using the “multiBatchNorm” function from the bachelor R package (v1.10.0). The per-gene variance of the log-expression profile was modeled using the “modelGeneVarByPoisson” function from the scran R package (v1.22.1), and the top 5,000 highly variable genes were selected. The mutual nearest neighbors (MNN) approach implemented by the “fastMNN” function from the bachelor R package was used to perform the batch correction. The first 50 dimensions of the MNN low-dimensional corrected coordinates for all cells were used as input to produce the t-stochastic neighbor embedding (t-SNE) projection and Uniform Manifold Approximation and Projection (UMAP) using the “runTSNE” and “runUMAP” functions from the scater R package (v1.22.0) respectively. Putative cell clusters were identified using the Leiden unsupervised community detection algorithm from the igraph R package (v1.3.0).

Marker genes and differential expression analysis

Cluster-specific markers were identified using the “findMarkers” function from the scran R package, which performs pairwise t-tests between clusters. Differential expression analysis between treatment and control (IgG) was performed on pseudo-bulk samples using the quasi-likelihood pipeline from the edgeR R package (v3.36.0). Genes with a false discovery rate below 5% were considered differentially expressed.

Pathway analysis

Pathway enrichment analysis of differentially expressed genes (DEGs) was performed using the R package enrichR (v3.0). Pathway results were ranked by the combined scores (calculated as log(p value) multiplied by z-score), and pathways with an adjusted p value<0.05 were considered significant. Further Gene Set Enrichment Analysis (GSEA) was done using the Metascape web tool12 and plotted as −log10P.

Subclustering and further analysis

Cells were divided into functional cell type groups and further subclustered using the Louvain method for unsupervised community detection.13 Trajectory inference was made using the Slingshot R package,14 and the association of our PARP14i signature with survival in Immune Checkpoint Blockade Therapy (ICBT) cohorts, as well as a comparison of high and low-expressed genes in specific clusters, was assessed using Tumor Immune Dysfunction and Exclusion (TIDE) web application.15 Sample RL1i projection on UMAP was done using the ProjecTILs R package.16

Data availability

Single-cell RNA sequencing (scRNAseq) data have been deposited in the ArrayExpress database at the European Molecular Biology Laboratory European Bioinformatics Institute (EMBL-EBI) under accession number E-MTAB-14471 (data set).17 Bulk RNA sequencing (RNAseq) was performed as previously described,5 with data available from the ArrayExpress database under accession numbers E-MTAB-12195 (data set)18 and E-MTAB-12872 (data set).19

ResultsA genetic signature reflecting the modulating effect of PARP14i on chronic IFNγ stimulation predicts response to α-PD-1 in clinical cohorts

We and others have previously implicated PARP14 in modulating type II interferon signaling.5 8 20 Following chronic stimulation of several mouse and human cancer cell lines with IFNγ and then subsequent treatment for 48 hours with a PARP14i, PARP14 PROTAC, or aqueous DMSO (drug solvent control), we assessed differences in gene expression by bulk RNAseq. Comparison of the top 50 DEGs within each cell model revealed a shared expression pattern, which we narrowed down to a PARP14i signature consisting of 30 genes (figure 1A, online supplemental figure 2A and online supplemental table 1). Pathway enrichment analysis performed by Metascape12 and ARCHS421 revealed that among this signature were genes strongly correlated with response to IFNγ signaling, with an emphasis on STAT/IRF function (figure 1B and online supplemental figure 2B). However, we found this signature is distinct from the genes contained within the IFNγ signaling Gene Ontology (GO) term, with less than half of the PARP14i signature displaying overlap. We then investigated whether this PARP14i signature could predict response to α-PD-1 therapy by employing the TIDE web application.15 High expression of genes within this signature was associated with substantial improvement in overall patient survival across several clinical cohorts receiving ICBT22–24 (figure 1C), underscoring the prognostic value of this signature and supporting the potential therapeutic benefit of PARP14 inhibition.

Figure 1

Identification of a genetic signature associated with PARP14 inhibition. YUMM2.1 and 501MEL cells were treated in vitro for two weeks with IFN-γ and then for 48 hours with either DMSO, 100 nM PARP14 inhibitor (RBN012759), or 500 nM PARP14 PROTAC (RBN012811); n=3. (A.) Bulk RNAseq analysis of a PARP14 inhibition signature identified across human and mouse cell lines. Heatmap of treatment versus control per cell line, colors representing row-scaled z-score expression high (red) to low (blue). Exp Factor indicates treatment with either DMSO (gray) or either PARP14i/PROTAC (green). (B.) Gene set enrichment analysis (GSEA) of 30 PARP14i signature genes identified by comparison of bulk RNAseq from human and mouse cell lines, chronically treated with IFN-γ followed by 48 hours of PARP14 inhibition (Metascape analysis, values presented as −log10 of p value). (C.) Kaplan-Meier (KM) curves showing patients with cancer with high expression of PARP14i signature genes having better overall survival (OS) rates in representative studies. Red indicates patients with high (top 50%) expression and blue indicates low (bottom 50%) expression. (D.) Bulk RNAseq analysis on naïve YUMM2.1 tumors treated with IgG or α-PD-1. α-PD-1-treated cohort was split into relapsed (REL) and responders (RES), compared with YUMM2.1 tumors pretreated with chronic IFN-γ stimulation. DMSO, Dimethyl sulfoxide; IFN-γ, interferon-gamma; IgG, immunoglobulin G; IRF1, interferon regulatory factor 1; PARP14, Poly ADP ribosyl polymerase 14; PD-1, programmed cell death 1; PD-L1, programmed cell death ligand 1; PROTAC, proteolysis targeting chimera; RNAseq, RNA sequencing; STAT1, signal transducer and activator of transcription 1.

Furthermore, we performed RNAseq analysis on subcutaneous tumors comprizing YUMM2.1 cells dissected from mice treated with either IgG or α-PD-1, dividing the latter into those either relapsing or responding to the treatment. Although the PARP14i signature was upregulated in α-PD-1-treated samples, it did not distinguish between treatment response and relapse (figure 1D). For further comparison, we also analyzed IgG-treated tumors arising from YUMM2.1 cells pretreated ex vivo with chronic IFNγ before implantation and found that the PARP14i signature was significantly downregulated (figure 1D), consistent with a tumor microenvironment (TME) that is T cell cold and inherently resistant to α-PD-1.5

Analysis of changes to the cellular composition of YUMM2.1 tumors relapsing following α-PD-1 therapy driven by subsequent PARP14 inhibitor treatment

To investigate further the ability of PARP14i to improve α-PD-1 response, we turned to our previously described mouse model of acquired resistance:5 YUMM2.1 cells were implanted subcutaneously, and on tumors reaching ∼80 mm3 in volume, two doses of α-PD-1 antibodies were given 3 days apart (figure 2A and online supplemental figure 3A). Tumor volumes were routinely recorded (online supplemental table 2), and at ∼100–150 mm3, animals were randomly allocated into the following retreatment groups: IgG2a (n=4); α-PD-1 (n=3); PARP14i alone (n=3); and a combination of α-PD-1 + PARP14 i (n=4). Animals were sacrificed following 7 days of retreatment, and tumors were subsequently processed for scRNAseq. Following tumor dissociation and enrichment for CD45+ cells, 5% from the CD45− cell fraction was added to ensure tumor cell representation in our analyses. Cells were fixed and frozen until all tumor samples were collected; then, samples were barcoded, multiplexed, and libraries prepared and sequenced.

Figure 2

Analysis of the cellular composition of YUMM2.1 tumors relapsing following α-PD-1 antibodies with and without PARP14 inhibitor treatment. (A) YUMM2.1 cells were implanted subcutaneously into the flank of female C57BL/6 mice. Once tumors reached ∼80 mm3 in volume, treatment of α-PD-1 antibodies was commenced for a total of two doses, 3 days apart. Tumor shrinkage and regrowth were monitored. When tumors reached a volume between ∼100–150 mm3, a second round of treatment commenced, lasting 7 days. Tumors were subsequently processed for single-cell RNAseq. Graphical summary of experimental timeline, showing tumor inoculation and treatment regime (top), and illustration of tumor growth (bottom). Created in BioRender. Leshem (2023) BioRender.com/r18g555. (B) Mutual nearest neighbor (MNN) corrected Uniform Manifold Approximation and Projection (UMAP) of 13 tumor samples (n=65,476 cells). Clustered using Leiden algorithm and colored by cell types. Cluster identifies below. (C) Cell number changes in each cluster by treatment: three mice IgG2a; three mice α-PD-1; three mice PARP14i (RBN012759) alone; and four mice a combination of α-PD-1 + PARP14 i (RBN012759). (D) Cell density projections by treatment. IgG, immunoglobulin G; NK, natural killer; PARP14, Poly ADP ribosyl polymerase 14; PD-1, programmed cell death 1; RNAseq, RNA sequencing.

In total, 14 samples of tumor immune cells were sequenced. Initial analysis of sequencing data identified one of the IgG samples (RL11i) to be anomalous, which on further review of its corresponding tumor growth curve was attributed to it being, in reality, a non-responder to the initial α-PD-1 treatment, and therefore, it was excluded from subsequent analyses (online supplemental figure 3A–C). Leiden unsupervised clustering across the remaining 13 treated samples identified 28 clusters based on gene expression patterns, with a diverse range of cell types represented across samples and treatment groups (figure 2B, online supplemental figure 3C,D and online supplemental tables 3 and 4). The following cell types were identified: Monocytes (Clusters 1–5); Macrophages (6–8); B cells (9–10); Plasma cells (11); Natural Killer (NK) cells (12); T cells (13–17); Dendritic cells (18–21); Granulocytes (22–25); Cancer-Associated Fibroblasts (CAFs) (26); and Tumor cells (27–28). To characterize the individual clusters within these groups, we manually validated cell type identity by comparing the expression of known markers across clusters (online supplemental tables 5–8).

Macrophages, Monocytes, and T cells were the most abundant leukocyte populations across all treatment groups, with distinct subpopulations displaying treatment-specific enrichment (online supplemental figure 3E). Notably, T cells were particularly abundant in combination-treated samples, with increased prevalence of clusters 14 and 16 and depletion of cluster 13, which was highly enriched in IgG-treated samples (figure 2C). Similarly, analysis by cell density plots revealed a shift from cluster 7 in macrophages, prominent in α-PD-1-treated samples, towards cluster 6 in the PARP14i and combination treatments (figure 2D). Furthermore, we observed IgG-specific enrichment and depletion of B cells and monocytes, respectively. Similarly, across our treatment cohort, few neutrophils and eosinophils were identified, and these could not be well separated, yet we found the sequenced tumor sample, which did not respond to α-PD-1 to be highly enriched in neutrophils (online supplemental figure 3B), consistent with a documented role of neutrophils in promoting rapid tumor growth.25

PARP14 inhibitor treatment increases M1 macrophage polarization to generate a memory-like state

Given the previous demonstration that PARP14i can suppress the reprogramming of tumor-associated macrophages (TAMs)26 that are widely considered to be protumorigenic27 and mimics the effects of α-PD-1 treatment on macrophages in renal tumor explants,7 we interrogated the effects of PARP14i on the macrophage populations in the TME within our specimens. By further subclustering clusters 1–8 from the initial analysis, we created a new MonoMac UMAP representation detailing a clear division of these subclusters into those that exhibited characteristics of classically activated (M1) macrophages (left) or alternatively activated (M2) macrophages (right) (figure 3A and online supplemental tables 6 and 9).

Figure 3

PARP14 inhibitor treatment increases M1 macrophage polarization to generate a memory-like state. (A) MNN corrected UMAP of monocytes and macrophages (clusters 1–8 of main data, n=34,705). Subclustered using the Louvain method of community ordering and colored by subclusters. (B) Cell number changes in each subcluster by treatment: IgG2a; α-PD-1; PARP14i (RBN012759); combination treatment α-PD-1 + PARP14i (RBN012759). (C) Cell density projections of both monocytes and macrophages (colloquially referred to as MonoMacs) by treatment. (D) Macrophage subcluster identity—subclusters 2, 4, 8 and 10. GSEA of upregulated (FDR=0.05, logFC>1, marked in red) and downregulated (FDR=0.05, logFC≤1, marked in blue) genes (Metascape analysis, values presented as −log10 of p value). (E) Cytokine and Chemokine gene expression presented as high (red) to low (blue) average expression and ratio of cells expressing (circle size), across Macrophage subclusters. (F) Gene expression presented as high (red) to low (blue) average expression and ratio of cells expressing (circle size), across Macrophage subclusters in different treatments. (G) Binarized heatmap of 10% random subsampling of each MonoMac subcluster by relative regulon activity (SCENIC analysis). FDR, false discovery rate; GSEA, Gene Set Enrichment Analysis; IgG, immunoglobulin G; IL 17, interleukin 17; IRF8, interferon regulatory factor 8; MNN, mutual nearest neighbors; PARP14, Poly ADP ribosyl polymerase 14; PD-1, programmed cell death 1; SCENIC, Single-Cell rEgulatory Network Inference and Clustering; STAT6, signal transducer and activator of transcription 6; UMAP, Uniform Manifold Approximation and Projection.

Analysis of cell type proportions within this subclustering revealed treatment-specific changes. Combination treatment notably shifted macrophage polarization from M2 towards M1 (figure 3B). Further inspection of gene enrichment analysis of DEGs, applied alongside cell density mapping, allowed us to robustly identify subclusters most affected by treatment (figure 3C,D). This demonstrated that subcluster 10—the most abundant cluster in the IgG-treated samples—is effectively depleted by combination treatment. This subcluster consists primarily of TAMs, which are known to facilitate tumor progression and metastasis, and as such, their reduction by combination treatment disrupts protumorous mechanisms.

These antitumoral effects are further supported by the expansion of specific macrophage subsets. We observed an increase in M1 macrophages (subcluster 2), signifying an enhanced proinflammatory response that bolsters antitumor immunity. This increase highlights the ability of combination treatment to promote a macrophage phenotype that is associated with effective immune surveillance and cytotoxic activity. Similarly, the increased abundance of a transitory M1/M2 macrophage state (subcluster 4) further illustrates the dynamic changes in macrophage populations, highlighting cellular plasticity modulated by combination treatment to influence the balance between inflammatory and immunoregulatory activities. Concurrently, macrophages in subcluster 8, most closely associated with an M2a macrophage identity, play a nuanced role in the TME under combination treatment. These macrophages exhibit high expression of Cxcl13 (figure 3E), which can mediate both protumorous and antitumoral effects. CXCL13 interaction with CXCR5 can recruit cytotoxic T cells to the TME, promoting tumoricidal factor secretion and supporting the formation of tertiary lymphoid structures.28 29 Conversely, this axis can also enhance IL-10 secretion by tumor cells, suppressing the proinflammatory functions of antigen-presenting cells (APCs) or recruiting immunosuppressive cells.30 In parallel, the expression of Ccl7 and Cxcl1 by M2a may also contribute to protumorigenic effects, such as promoting angiogenesis and invasion,31 32 reflecting the dual functionality of this macrophage subset in the TME.

The overall increase in M1 macrophages following combination treatment is linked to elevated expression of proinflammatory genes (eg, Tnf and Stat3) (figure 3F), particularly in subclusters 2 and 4, while the decreased expression of Stat6 and Jun correlates with their reduced transcriptional activity and cellular activation. These changes suggest subclusters 2 and 4 may represent memory-like cells with the potential for reactivation.33 In contrast, PARP14i treatment results in a marked downregulation of Nfatc1, Nfatc3, and Smad4 across subclusters, reducing Tgfb1 signaling in M2 macrophage clusters, as evidenced by lower Sox2 and Sox4 expression (figure 3F).34 The depletion of key drivers of M2 macrophage polarization such as Tgm2, Vegfa, and Cebpb35 in clusters 8–11 by combination treatment aligns with an overall reduction in M2 macrophage populations. In the context of subcluster 8, however, which is enriched by combination therapy, the reduced levels of anti-inflammatory gene expression by these cells (figure 3F) positions them as a transitional macrophage state with reduced features of classical M2 macrophages, reinforcing a combination treatment-driven transition towards an anti-tumoral TME.

To further explore the molecular underpinnings of macrophage polarization, we applied Single-Cell rEgulatory Network Inference and Clustering (SCENIC) analysis36 to evaluate transcriptional regulatory programs within these subclusters (figure 3G). Here, we identified regulons that effectively delineated M1 (Runx3 and Cux1) and M2 (Irf4) macrophages. Additionally, six regulons associated with M2 macrophages (Nfia, Nr1d1, Klf2, Klf3, Etv1, Dbp) were most strongly associated with subcluster 10. In the context of combination treatment, the depletion of subcluster 10 resulted in a significant reduction in the activity of these regulons. However, this reduction may be offset by the enrichment of subcluster 8, indicative of a directional shift in this treatment towards an M1-like phenotype. These findings suggest that the overall balance between M1 and M2 macrophages is shifting towards a proinflammatory, antitumor state in combination treatment. Notably, subclusters 2 and 4, despite their relatively low regulon activity beyond the M1-defining features, may represent more memory-like macrophage states.

PARP14 inhibitor treatment in combination with α-PD-1 therapy drives activated cytotoxic T cells and limits their transition to dysfunctional states

PARP14 (also known as ARTD8) has been shown to interact with STAT6, a transcription factor activated by IL-4, to promote Th2 cell differentiation,37 while Hoon Cho et al demonstrated that PARP14 is involved in the regulation of antibody responses, particularly IgA production, through promoting Th17 differentiation.38 Both Th2 and Th17 cells are implicated in immune dysfunction in tumors.39 Moreover, our previous work demonstrated that chronic PARP14i treatment of stimulated T cells resulted in increased percentages of IFNγ-producing and TNFα-producing CD4 and CD8 T cells, accompanied by decreases in TGFβ and IL-10,5 consistent with a shift toward a type 1 proinflammatory T cell response in the TME.

To further characterize specific T cell identities and assess their functionality, we conducted subclustering of clusters 12–17 from the initial clustering (figure 4A). By combining NK and T cells for further subclustering, we can be confident that our new T cell clusters are free from any contaminating readouts due to the presence of NK T cells. This subclustering can be split broadly into NK and NKT cells (1–5) and T cells (6–19), with key subpopulations including naïve T cells (7 and 19), regulatory T cells (10 and 11), and cytotoxic T cells (13–18). Comparing the proportions of each T cell subcluster revealed treatment-specific changes in cellular identity and function (figure 4B and online supplemental table 10).

Figure 4

PARP14 inhibitor treatment in combination with α-PD-1 therapy drives activated cytotoxic T cells but limits transition to dysfunctional states. (A) MNN corrected UMAP of T cell and natural killer (NK) cells (clusters 12–17 main data, n=18,298). Subclustered using the Louvain method of community ordering and colored by subclusters. (B) Cell number changes in each T and NK cell subcluster by treatment: IgG2a; α-PD-1; PARP14i (RBN012759); combination treatment α-PD-1 + PARP14i (RBN012759). (C) Cell density projections of T and NK cells by treatment. (D) Cytokine and Chemokine gene expression presented as high (red) to low (blue) average expression and the ratio of cells expressing (circle size), across T and NK cell subclusters. Subclusters of interest highlighted on y-axis. (E) T cell subcluster identity—clusters 13–16. GSEA relevant terms of upregulated (red, FDR=0.05, top 25 genes) and downregulated (blue, FDR=0.05, bottom 25 genes) genes (Metascape analysis, values presented as −log10 of P value). (F) Slingshot trajectory inference analysis of CD8+ T cell subclusters 13–18 (n=6,194) with trajectory direction indicated by arrows (black). (G) KM curves showing cancer patients with high expression of genes from T cell subclusters 13–16 and their effect on OS rates in bladder cancer cohort. Red indicates patients with high (top 50%) expression and blue indicates low (bottom 50%) expression. (H) Binarized heatmap of each T cell subcluster by relative regulon activity (SCENIC analysis). FDR, false discovery rate; GSEA, Gene Set Enrichment Analysis; IgG, immunoglobulin G; IRF4, interferon regulatory factor 4; KM, Kaplan-Meier; MNN, mutual nearest neighbors; OS, overall survival; PARP14, Poly ADP ribosyl polymerase 14; PD-1, programmed cell death 1; PD-L1, programmed cell death ligand 1; SCENIC, Single-Cell rEgulatory Network Inference and Clustering; Treg, regulatory T cell; UMAP, Uniform Manifold Approximation and Projection.

Cell density mapping (figure 4C) highlighted a treatment-dependent decrease in naïve T cells (subcluster 7), which were highly abundant in IgG-treated samples, suggesting their transition into differentiated or activated states on treatment. In accordance, in IgG-treated tumors, naïve T cells (subcluster 7) exhibited high expression of Ccr7 and Ccr9, which likely mediated their recruitment to the TME (figure 4D). CCR7 facilitates migration toward lymphoid tissues via CCL19, a chemokine produced by exhausted CD8+ T cells (subcluster 13) and M2a macrophages (MonoMac subcluster 8).40 This interaction may reflect a futile attempt to restock the TME with tumor-inexperienced T cells. However, these naïve T cells exhibited increased Sell expression compared with other T cell subclusters, indicative of early activation (online supplemental figure 4A) and that IgG-treated tumors suffer from dysfunctional T cell priming by APCs.

Corroborating our hypothesis that the reduction in naïve T cells is due to their differentiation into more active T cell states, we identified a significant enrichment of effector T cells (subclusters 13–16) in the combination treatment group. These cells are recruited and retained in the TME through interactions with proinflammatory chemokine ligands, including CCL3, CCL4, and CCL5 (figure 4D).41 This recruitment is further reinforced by the CXCR5-CXCL13 axis, where CXCL13, predominantly produced by M2a macrophages in our data as previously described (MonoMac subcluster 8), drives cytotoxic T cell recruitment to the TME and aids the formation of tertiary lymphoid structures.28 29 These structures serve as immune activation hubs, enabling CXCR5-expressing effector T cells to undergo further activation and proliferation, enhancing their cytotoxic function.

Subclusters 13–16 display distinct functional profiles associated with T cell exhaustion, quiescence, hyperactivation, and T cell response, respectively (figure 4E and online supplemental figure 4B,C). In the context of antitumor activity, T cell subclusters 13–16 are closely interconnected in driving the T cell effector response. Our slingshot trajectory analysis (figure 4F) indicates that by entering a state of dormancy (subcluster 14), hyperactive T cells (subcluster 15) sidestep terminal exhaustion (subcluster 13) in the presence of chronic antigen stimulation. Thereafter, quiescent cells can either undergo reactivation (subcluster 16) or persist as a reservoir of memory T cells with retained effector function, poised for this reactivation.

Kaplan-Meier survival plots were constructed using the TIDE web application for gene expression signatures arising from each T cell subcluster from subclusters 13–16 (online supplemental table 11), revealing that high expression of the gene signature for subcluster 16 (reactivated T cells) was associated with significantly improved survival of patients receiving ICBT for either melanoma or bladder cancer (figure 4G and online supplemental figure 4D). Furthermore, we observed this was also true, although to a lesser degree, for the gene signatures arising from subclusters 14 (quiescent T cells) and 15 (hyperactive T cells). On the other hand, the gene signature of T cell subcluster 13 (exhausted T cells) did not associate with positive patient outcomes, in keeping with T cell exhaustion constituting a profound loss of T cell functionality.42 This suggests that cells within subclusters 14–16 are critical in mounting an effective antitumor response through mechanisms that sustain effector function without succumbing to terminal exhaustion.

We next employed SCENIC analysis to highlight transcriptional regulatory systems that are most pertinent to the functional and prognostic characteristics of these T cell populations and revealed distinct regulatory patterns across the T cell subclusters, highlighting the differential activity of several key transcription factors (figure 4H and online supplemental figure 4E). Jdp2, an AP-1 repressor, exhibited differential regulon activity across the T cell subclusters, being elevated in subcluster 13 compared with subclusters 14–16. Higher Jdp2 expression contributes to impaired cellular functionality by T cell exhaustion and anergy.43 44 Therefore, lower activity of the Jdp2 regulon in subclusters 14 and 15 is consistent with a more functional and active T cell state and, considering the enrichment of these subclusters with combination treatment supports the notion that decreased Jdp2 activity facilitates better immune responses. Furthermore, elevated activity of the Mxd1 regulon in subcluster 13 suggests this, too, has a crucial role in maintaining the exhausted state characteristic of this subcluster by repressing cellular growth, contributing to its adverse clinical outcomes.