Prior anti-CTLA-4 therapy impacts molecular characteristics associated with anti-PD-1 response in advanced melanoma

Immune checkpoint inhibitors (ICIs), including CTLA-4- and PD-1-blocking antibodies, can have profound effects on tumor immune cell infiltration that have not been consistent in biopsy series reported to date. Here, we analyze seven molecular datasets of samples from patients with advanced melanoma (N = 514) treated with ICI agents to investigate clinical, genomic, and transcriptomic features of anti-PD-1 response in cutaneous melanoma. We find that prior anti-CTLA-4 therapy is associated with differences in genomic, individual gene, and gene signatures in anti-PD-1 responders. Anti-CTLA-4-experienced melanoma tumors that respond to PD-1 blockade exhibit increased tumor mutational burden, inflammatory signatures, and altered cell cycle processes compared with anti-CTLA-4-naive tumors or anti-CTLA-4-experienced, anti-PD-1-nonresponsive melanoma tumors. We report a harmonized, aggregate resource and suggest that prior CTLA-4 blockade therapy is associated with marked differences in the tumor microenvironment that impact the predictive features of PD-1 blockade therapy response.

Immune checkpoint inhibitors (ICIs) have revolutionized melanoma treatment. The anti-CTLA-4 agent ipilimumab was approved in 2011, followed by approval of the anti-PD-1 agents pembrolizumab and nivolumab in 2014. 1,2 With 5-year overall survival rates up to 52% in melanoma with combinatorial ICI approaches, 3 there is substantial effort to identify predictive biomarkers for response and resistance. Correlative studies have demonstrated patient stratification because of genomic and transcriptomic features; however, few discoveries, despite vali-dation within and across cohorts, have been successfully prospectively applied to clinical trials.
Anti-PD-1 therapies have been approved for patients with tumors with mismatch repair deficiencies, leading to increased tumor mutational burden (TMB), 4 or for solid tumors with high TMB (over 10 mutations/Megabase [mut/Mb]), following progression on prior treatments and without additional treatment options. 5 However, TMB and its inferred neoantigen burden have shown inconsistent concordance with response across and within tumor types and treatment settings. [6][7][8] Bulk gene expression and spatial profiling for tumor immune infiltration or exclusion (e.g., TIDE, IMPRES), interferon gamma (IFNG) pathway activation (e.g. GEP), antigen presentation machinery, and antitumor immune cell populations have also been used to positively identify patients who will benefit from PD-1 blockade, [9][10][11][12][13][14][15] and combination of these patterns with high TMB has been shown to be even more effective at predicting anti-PD-1 response. 7,12 As more correlative datasets are generated, there have been challenges in discovering and validating ICI biomarkers. Recent studies have attempted to overcome limitations in statistical power by aggregating molecular datasets. 7,12,16,17 However, in doing so, it is important to note the level of heterogeneity in the dataset and whether the annotation of technical, biological, and clinical variables is sufficient to delineate appropriate biomarkers and clinical contexts for application. Technical batch effects may include sample collection and preservation, nucleic acid extraction, sequencing approach, and methods for analysis. Additional biological challenges in these studies include the fundamental differences in the mechanism of each ICI, intrapatient or disease heterogeneity in genetic drivers, and immunogenicity within and across cancer types. Within melanoma, clinical trials are often inclusive of multiple melanoma subtypes, which have known differences in genomic profiles and observed differences in response to ICIs. 8, 18 Acral and mucosal melanomas have demonstrated reduced response and decreased survival compared with cutaneous and occult (unknown) melanomas and have reduced TMB, reduced UV-induced mutagenesis, fewer BRAF and NRAS hotspot mutations, and more frequent KIT mutations. 8, 18,19 Here, we present a harmonized dataset of whole exome sequencing (WES) and RNA sequencing (RNA-seq) samples derived from melanoma tumor biopsies collected from 514 subjects treated with ICIs enrolled across seven clinical cohorts, two of which have not been reported previously that are from the CheckMate 064 and 067 clinical trials. 20,21 We demonstrate that efforts to aggregate and harmonize molecular and clinical annotation of datasets facilitates statistical detection of genomic and expression-based correlates of response. Furthermore, exploration of clinical factors, such as melanoma subtype and prior immunotherapy treatment, reveals that previously reported predictors of anti-PD-1 response (e.g., TMB, types of immune infiltration) vary based on prior anti-CTLA-4 exposure. We summarize the collective immune infiltrate of biopsies using an aggregate metric, tumor immunogenicity associated with response to anti-PD-1 (TIARA-PD-1), to further define tumor-intrinsic and microenvironment expression patterns related to the overall immune status of each biopsy. Our work provides a single, cohesive genomic reference dataset for uncovering molecular characteristics associated with ICI response and resistance in melanoma, and it may have future utility for exploring molecular and clinical stratification of this disease.

Cohort summary
We aggregated seven clinical cohorts of patients with melanoma treated with ICI regimens, comprising 514 patients with either WES (N = 339), RNA-seq (N = 403), or matched WES and RNA-seq (N = 258; Figure 1A) derived from melanoma tumor biopsies. 6,8,[20][21][22][23][24] While some of these datasets have been published previously, 22% of the patients enrolled in these trials (N = 115 patients) 20,21 with correlative WES (n = 184 tumors) or RNA-seq (n = 95; n = 58 paired) have never been interrogated, providing an additional rich source of data. The following demographics were assembled for all patients and corresponding (B) Alluvial plot depicting the demographics of patients with melanoma tumor samples in the final dataset (x axis; cohort, subtype, prior ICI therapy, treatment regimen, and RECIST response). Each individual (alluvium) is colored by whether the corresponding sample was included in analysis of cutaneous melanoma tumors treated with anti-PD-1. See Table 1 and Figure S1. tumor specimens: RECIST 1.1 response to therapy, melanoma subtype, sex, age, treatment regimen, prior CTLA-4 blockade therapy status, originating study cohort, and time point with respect to therapy ( Figure 1B; Table 1). Treatment regimens included anti-PD-1 monotherapy, anti-CTLA-4 monotherapy, concurrent combination of these two therapies, or sequential treatment with these two therapies (i.e., anti-PD-1 to anti-CTLA4 or anti-CTLA4 to anti-PD-1).
Somatic variant calling and gene expression quantification were harmonized across all cohorts ( Figure S1), and tumornormal WES pairs were evaluated for a range of quality control metrics (Figures S1 and S2; STAR Methods). After samples were excluded because of low sequencing coverage, low estimated tumor purity, or the presence of unmatched sample contamination, the final cohort presented in the following analyses included 427 patients, 284 tumor biopsies with WES, and 442 biopsies with RNA-seq (188 biopsies had matched WES and RNA-seq data). With these efforts, we aimed to establish a publicly available resource summarizing these data and to utilize the scale of this dataset to specifically evaluate genomic and transcriptomic correlates of anti-PD-1 response in cutaneous melanoma.
We were specifically interested in the association of TMB and response to anti-PD-1 therapy in patients with cutaneous melanoma, leveraging the statistical power of our cohort. When TMB was compared across biopsies of cutaneous melanoma origin treated with anti-PD-1 that responded (RECIST complete or partial response [CR/PR]) or progressed (RECIST progressive disease [PD]), TMB was significantly greater in biopsies of patients with anti-PD-1 clinical response compared with no response (Wilcoxon test, p = 0.024). This association only remained significant, however, in tumors from patients that were previously treated with anti-CTLA-4 therapy (Wilcoxon test, p = 0.0062, FDR = 0.037) and not in anti-CTLA-4-naive tumors (Wilcoxon test, p = 0.515, FDR = 0.51) ( Figure 2B), suggesting an interaction effect between these two variables. This was consistent with what has been reported previously, 8 while other subsets were not of sufficient size to identify these differences ( Figure S2).
Single-nucleotide variants in PIK3C2G are positively associated with clinical benefits to checkpoint inhibitors While studies have attempted to discover associations between response to immune checkpoint blockade and single gene mutational statuses, success has been partially limited by the larger patient numbers required for statistical power following multiple hypothesis correction. 27 Tumor sampling and purity, prior treatment, and melanoma subtype, which are also associated with prognosis following ICI and the genomic landscape, 8 add complexity to this problem. We explored whether our aggregation and harmonization may overcome these challenges by applying logistic regression to genes and mutation types in cutaneous melanoma tumors treated with anti-PD-1, comparing results from biopsies of patients with response (CR/PR) with no response (PD). Gene-mutation combinations were filtered to those that occurred at least 10 times (of 151 [6.6%] tumors evaluated), resulting in 7,579 unique gene-mutation type pairs evaluated. There were two nonsilent alterations that met the statistical thresholds following multiple testing correction (p < 0.05, FDR < 0.2), enrichment of PIK3C2G missense mutations (p = 7.5eÀ5, FDR = 0.15), and DNAH5 splice region variants (p = 6.2eÀ3, FDR = 0.18) in cutaneous tumors from patients with response to anti-PD-1 ( Figure 2C). When we repeated this analysis across all melanoma subtypes and ICI regimens, PIK3C2G missense mutations remained associated with anti-PD-1 response across all melanoma subtypes (p = 1.3eÀ5, FDR = 0.033) as well as response to any ICI regimen across all subtypes (p = 3.2eÀ6, FDR = 0.011).
The most common alterations in PIK3C2G included the recurrent missense mutations R779C and E1272K ( Figure 2D). Of the 49 tumors with PIK3C2G mutations, 44 had matched RNA-seq data; however, when mutations were queried in RNA-seq, the gene and the mutation were infrequently expressed. The expression of PIK3C2G was low across all available melanoma tumors with RNA-seq data (0-14.68 fragments per kilobase transcript per million mapped reads [FPKM]; median 0). Notably, the link between PIK3C2G mutations and anti-PD-1 response has been demonstrated previously but not specifically reported because it was not significant following multiple hypothesis correction. 8 The inclusion of additional harmonized cohorts (CheckMate 038, 064, and 067) ( Figure 2E) provided additional key support for this finding.
Prior CTLA-4 blockade differentially stratifies the expression landscape of melanoma tumors with and without response to anti-PD-1 Previous studies have described remodeling of the tumor microenvironment (TME) by anti-CTLA-4 blockade. 28,29 Because of the differences observed in TMB with respect to prior anti-CTLA-4 therapy and anti-PD-1 response, we performed supervised differential expression analysis between baseline cutaneous melanoma biopsies from anti-CTLA-4-experienced patients (n = 39; n = 17 CR/PR, n = 22 PD) and anti-CTLA-4-naive patients (n = 92; n = 44 CR/PR, n = 48 PD) treated with anti-PD-1 monotherapy ( Figure 3A). When patients were stratified by their response to anti-PD-1, we found more significant differences across biopsies of patients with response to anti-PD-1 (CR/PR; n = 1,025 significant genes) than across biopsies of patients  Genes that were differentially expressed between anti-CTLA-4-experienced and anti-CTLA-4-naive tumors were mostly unique to either response or no response to anti-PD-1 therapy, with 63.5% and 48.7% of genes specific to either subset, respectively ( Figure 3D). The lack of overlap in these results suggested that there were biological differences between melanomas responsive and nonresponsive to PD-1 blockade, based on their prior treatment with anti-CTLA-4 ( Figure 3E). Only 11 genes were shared by both analyses; 9 of 11 genes had shared direction of effect.
Given the divergence in gene expression patterns between anti-CTLA-4-experienced and anti-CTLA-4-naive tumors within each anti-PD-1 response group, we utilized two methods for gene set enrichment analysis (GSEA) to more broadly characterize molecular processes across these comparisons. First, GSEA was performed within each response group using only significantly different genes (FDR < 0.10). Consistent with our results at the gene level, there were stronger differences at the gene set level when comparing anti-CTLA-4 experience in anti-PD-1-responsive tumors using Gene Ontology (GO): biological process (BP) gene sets (Figure 3F). In anti-PD-1-responsive tumors, those with anti-CTLA-4 experience were significantly enriched for processes related to mitogen-activated protein kinase (MAPK) signaling, hypoxia response, and inflammatory signaling pathways (e.g., immunoglobulin G [Ig] production, IFN signaling, adaptive immune response), while anti-CTLA-4-naive tumors were enriched for cell-cycle pathways. The same analysis revealed very few significantly enriched pathways in biopsies of patients with no response to anti-PD-1, with two exceptions: the adaptive immune response gene set was enriched in anti-CTLA-4-experienced tumors, while the innate immune system gene set was enriched in anti-CTLA-4-naive tumors.
In our second GSEA method, we performed a rank-based GSEA in Hallmark gene pathways 30 across all genes from the prior differential gene expression analysis. Within the anti-PD-1-responsive subset, anti-CTLA-4-experienced tumors were enriched for gene sets associated with IFNG response and genes downregulated by KRAS activation, while anti-CTLA-4-naive tumors were enriched for cell cycle, E2F, G2M checkpoint, and mTORC1 signaling pathways ( Figures 3G and 3H). Nearly all of the significant gene sets in anti-PD-1-responsive tumors showed  (legend continued on next page) opposing, but not significant, trends in anti-PD-1-nonresponsive tumors (e.g., gene sets enriched in anti-CTLA-4-naive tumors in the anti-PD-1-responsive subset were enriched in anti-CTLA-4experienced tumors in the anti-PD-1-nonresponsive subset) ( Figures 3G and 3H). The only two shared significant processes across anti-PD-1 responding and non-responding tumors were enrichment of fatty acid metabolism and myogenesis in anti-CTLA-4-experienced tumors.
Collectively, these two approaches revealed interactions between anti-CTLA-4 experience and anti-PD-1 response, particularly in responsive tumors. We highlight that the majority of pathways from GSEA would not have achieved statistical significance without the stratification of anti-PD-1 baseline tumors by clinical response, supporting continued stratification of patients by anti-PD-1 response and prior anti-CTLA-4 therapy for the subsequent analysis.
The strength of inflammatory gene expression patterns in anti-PD-1 responding biopsies is associated with prior anti-CTLA-4 experience To more specifically explore gene expression patterns associated with anti-PD-1 response, we stratified cutaneous melanoma tumors with RNA-seq data (n = 131) into four subsets based on anti-PD-1 response (CR/PR, PD) and prior anti-CTLA-4 experience (experienced, naive). Unsupervised principal-component analysis of the top 15% most variable genes did not reveal separation across clinical demographics ( Figures 4A, S3A, and S3B). We utilized four models for supervised differential expression analysis to compare biopsies of anti-PD-1 response with no response: comparison of all samples, comparison of all samples controlling for anti-CTLA-4 experience, comparison of response within anti-CTLA-4-experienced tumors, and comparison of response within anti-CTLA-4naive tumors ( Figures 4B-4D). By comparing the results of these approaches, we corroborated our previous results indicating interactions between anti-CTLA-4 experience and anti-PD-1 response ( Figure 4E). This was supported by analysis of each patient subset, where there were no significantly different genes comparing anti-PD-1 response in anti-CTLA-4-experienced tumors.
We explored whether differences were driven by one or a subset of clinical cohorts, and we found that the differential expression analysis within subgroups was poorly powered ( Figure S3C). When the ranks of significant genes were compared between the analysis of the aggregate dataset with each individual cohort, positive correlations were significant but weak (Pearson correla-  Figure S3D). Notably, the Liu et al. 8 cohort did not have the strongest rank-based correlation despite having the largest sample contribution. This sensitivity analysis underscores the importance of statistical power provided by our data harmonization and aggregation efforts to understand associations between molecular and clinical features.
We further explored these gene patterns by performing GSEA across all genes of curated Hallmark and KEGG pathways, revealing striking differences between anti-CTLA-4-experienced and anti-CTLA-4-naive tumors when comparing biopsies of patients with and without response to anti-PD-1 (Figures 4F and  4G). Within anti-CTLA-4-experienced tumors, anti-PD-1 responding tumors were enriched for genes regulated by KRAS activation, while anti-PD-1 non-responding tumors were enriched for cell cycle and MYC target genes. However, these gene sets displayed the reverse patterns in anti-CTLA-4-naive tumors ( Figure 4G). Gene sets related to inflammation, cytotoxicity, and immune cell activation were positively enriched in anti-PD-1-responding tumors, but enrichment was significant only in anti-CTLA-4-experienced patients. The expression of cytokines and immunoglobulin regions was also variable, with the highest expression in anti-CTLA-4-experienced tumors of patients who responded to anti-PD-1, followed by anti-CTLA-4naive without a response, with the lowest expression in biopsies of patients with no response to anti-PD-1 ( Figure 4H). Consistent with prior studies, these results suggest that biopsies of patients with response to anti-PD-1 have increased immune infiltration and activity compared with non-responding biopsies, but these features are more strongly enriched in biopsies of patients who are anti-CTLA-4 experienced and responded to anti-PD-1.
Anti-CTLA-4 experience is associated with globally increased adaptive immune infiltration signatures While these analyses revealed differences in expression or enrichment for inflammatory processes associated with anti-CTLA-4 experience and anti-PD-1 response, we were interested in which immune cell types may be contributing to these signals. Previous studies have described the association between tertiary lymphoid structures (TLSs) and infiltration of antitumor CD8 T cells and anti-PD-1 response 9 as well as the mechanism by which anti-CTLA-4 mediates T cell trafficking from the lymph node to the tumor site. 28,29 To interrogate immune cell types in the TME, we computed a single-sample immune signature score using gene signatures from previously published and curated gene sets to estimate immune cell type enrichment and activity from bulk tumor RNA-seq data (STAR Methods). Concordant with the gene expression analysis, anti-CTLA-4-experienced, anti-PD-1-responsive tumors had the greatest average enrichment across all immune cell signatures, with the strongest enrichment for lymphocyte, conventional type 1 dendritic cells (cDC1), and TLS gene signatures ( Figure 5A). Cell signature estimation using EPIC demonstrated similar trends with anti-CTLA-4 experience. 31 The EPIC endothelial signature also suggests increased blood vessels in the CR/PR anti-CTLA-4-experienced group compared with PD groups (Figures S4A and S4B).
To determine which immune cell types may be cooperative within the TME, we normalized the raw cell type signatures for each immune cell type within each sample and summarized the average score for each cell type within each clinically defined subset. This revealed that anti-CTLA-4-experienced,   Figures 5A and 5B). Conversely, anti-PD-1 non-responding biopsies had greater enrichment of Treg cell, neutrophil, and MDSC signatures. While anti-CTLA-4-experienced tumors appear to have stronger overall enrichment of immune cell signatures, the anti-PD-1-non-responding biopsies also showed reduced lymphocyte infiltration and TLS signatures along with increased NK cell signatures. When we compared the average enrichment score of each cell-type between anti-PD-1 responding and non-responding biopsies, there was an overall inverse association in cell type signatures across all cell-types ( Figure 5C). Importantly, this difference across anti-PD-1 response groups was more substantial in anti-CTLA-4-experienced tumors.
Tumor infiltration is differentially associated with MYC and E2F targets in biopsies that responded to anti-PD-1 after anti-CTLA-4 treatment Gene expression patterns associated with immune infiltration, inflammatory signaling, and antitumor immune cytotoxicity have been shown to positively predict response to anti-PD-1 from baseline tumor specimens. Several methods have been established to aggregate expression metrics across gene sets from bulk RNA data, leading to successful stratification of anti-PD-1 or ICI responders. 10,13-15 However, we aimed to understand how tumor-intrinsic patterns may coordinate with immune cell types responsible for anti-PD-1 response. To better stratify tumors based on the immune infiltrate, we generated an interpretable, aggregate gene expression score, using immune cell type signatures that were increased or decreased in biopsies of patients with response to anti-PD-1, to summarize the TIARA-PD-1 in melanoma (STAR Methods).
Given that the features defining TIARA-PD-1 were computed based on anti-PD-1 response in our aggregate cohort of baseline, cutaneous melanoma tumors, it was unsurprising that the TIARA-PD-1 score was significantly increased among biopsies of patients with CR/PR to anti-PD-1 compared with PD (Figure 5D). Consistent with our previous analysis, this difference was greater in magnitude when including prior anti-CTLA-4 treatment as a covariate in biopsies of patients with anti-PD-1 response but not among nonresponsive tumors. We also computed TIDE, estimating T cell dysfunction and exclusion, and IMPRES, a predictor of ICI response in melanoma that includes 15 pairwise transcriptomics relations between immune checkpoint genes. 13,15 TIDE and IMPRES support interactions with prior anti-CTLA-4 experience and response to anti-PD1 ( Figure S4C). However, while TIARA-PD1 variance is associated with CR/PR in the naive and anti-CTLA-4 experienced groups, TIDE was only able to stratify anti-PD-1 response in the anti-CTLA-4-experienced setting ( Figure S4D). Importantly, compu-tation of TIDE scores required manual supervision of whether the patients had experienced prior immunotherapy to identify these interactions (STAR Methods).
We sought to validate whether TIARA-PD-1 was associated with clinical response in melanoma tumors in other ICI contexts by quantifying immune cell type signatures and TIARA-PD-1 across baseline and on-treatment RNA-seq samples from tumors treated with other ICI regimens associated with this study (Figures S5A and S5B). Consistent with our anti-PD-1 results in baseline cutaneous melanoma, clinical response (CR/PR) to ICI in other treatment regimens tended to have larger TIARA-PD-1 scores ( Figure S5C). Additionally, median TIARA-PD-1 was higher in on-treatment samples compared with baseline in the anti-PD-1 and sequential anti-CTLA-4-to-anti-PD-1 treatment regimens, suggesting increased global immune infiltration ontherapy with ICI. We also performed an independent validation of TIARA-PD-1 in a small clinical trial in which patients with advanced melanoma were treated with either anti-CTLA-4 monotherapy or anti-CTLA-4 plus anti-PD-1 combination therapy. 32 Tumors with either PR or stable disease (SD) in this cohort did show a trend of increased TIARA-PD-1 over the course of therapy, comparing lesion-paired on-treatment samples (n = 10) with baseline samples (n = 10) (Figures S5D and S5E). When we evaluated TIARA-PD-1 in melanoma subtypes other than cutaneous melanoma, TIARA-PD-1 was associated with anti-PD-1 response in mucosal melanoma but not in tumors with acral, uveal, or unknown site of origin ( Figure S5F). Notably, these comparisons were less powered than the analysis of cutaneous melanomas.
To better understand the associations between the TME and potential tumor-intrinsic processes, we performed differential expression analysis in baseline cutaneous melanoma tumors treated with anti-PD-1 based upon TIARA-PD-1 score as an outcome variable. Gene expression analysis identified 3,961 significantly different genes (FDR % 0.1; Figure 5E). Genes contributing to the immune signatures used to compute TIARA-PD-1 were excluded from this analysis. Unsurprisingly, increased TIARA-PD-1 was positively associated with LCK, immunoglobulins, XCL1, and other cytokines ( Figure 5E), and GSEA on ranked genes revealed enrichment for Hallmark gene sets associated with T cell function, IFN signaling, and tumor necrosis factor alpha (TNF-a) signaling via nuclear factor kB (NF-kB) across all clinical subsets ( Figure 5F; Table S5). Conversely, TIARA-PD-1 was negatively associated with myosin-and actinfamily genes, with enrichment of genes associated with myogenesis, epithelial-to-mesenchymal transition, glycolysis, MYC targets, and hypoxia ( Figures 5E and 5F).
TIARA-PD-1 was informed by the cell types enriched in anti-PD-1 responding biopsies, and our prior analysis showed that this signal was greater in anti-CTLA-4-experienced tumors. Differential gene expression analysis of each clinical subset showed the greatest number of significant genes associated (E) Differential expression analysis as a function of TIARA-PD-1. The x axis is the expression difference associated with a unit change in TIARA-PD-1. (F) NES of rank-based GSEA for statistically significant pathways curated from Hallmark and NABA (Table S5) Figures S4-S6 and Tables S4-S7. with TIARA-PD-1 among anti-CTLA-4-naive responders, followed by anti-PD-1 non-responders, anti-CTLA-4-experienced non-responders, and then anti-CTLA-4-experienced responders ( Figure 5G). The low number of significant genes associated with TIARA-PD1 in the anti-CTLA-4 experienced subsets may be a result of the smaller sample size, as seen in our previous transcriptomics analysis. When GSEA was performed within each clinical subset, many gene sets were consistently enriched in the same direction, particularly those that were negatively associated with TIARA-PD-1 ( Figure 5H; Table S6). Gene sets that had opposing directions in enrichment scores were only significant in one or two clinical subsets. In particular, MYC and E2F targets showed diverging associations with TIARA-PD-1 in anti-PD-1-responsive tumors, with positive enrichment in anti-CTLA-4-experienced tumors and negative enrichment in anti-CTLA-4-naive tumors. MYC targets were also negatively enriched in anti-CTLA-4-experienced, anti-PD-1-non-responding tumors ( Figure 5H). Collectively, these results suggest differential roles in these processes in either coordinating with or facilitating antitumor immune activity responsible for anti-PD-1 response.

Mutation subtype may be associated with immune cell infiltration and anti-PD1 response in an anti-CTLA-4experience-dependent manner
We sought to further understand whether the patterns with respect to prior anti-CTLA-4 and anti-PD1 response observed in baseline cutaneous melanoma samples were associated with known canonical genetic drivers of melanoma (BRAF hotspot, NRAS hotspot, NF1 mutant, or triple wild type [TWT]). Driver mutations detected in WES data were annotated based on detection of hotspot mutations in BRAF (V600) or NRAS (G12, G13, Q61), mutations in NF1, or the lack of these mutations (TWT). Clinical benefit to anti-PD1 or prior anti-CTLA4 experience was not associated with any genetic subtype (two-sided Fisher's exact test, p = 0.53). We performed differential expression analysis comparing anti-PD-1 response (CR/PR versus PD) for anti-CTLA-4-naive and anti-CTLA-4-experienced tumors within each genetic subtype. There were no significant genes identified in anti-CTLA-4-naive tumors. In the anti-CTLA-4-experienced group, IGHM was the only significant gene positively associated with CR/PR in tumors with BRAF hotspot mutations (log2 fold change = 3.56, FDR = 0.099; Table S7). Notably, sample sizes were small after stratification by anti-PD-1 response, anti-CTLA-4 experience, and mutation subtype.
Regardless of genetic subtype, anti-CTLA-4-experienced patients who responded to anti-PD1 tended to have greater global immune infiltration than those without benefit ( Figure S6A). There was no significant association between immune cell signatures across mutant subtypes or within clinical groups, stratified by mutant subtypes. TIARA-PD1 analysis of immune signatures was consistent across genetic subtypes, with the highest TIARA-PD1 in anti-CTLA-4-experienced, anti-PD-1 responders ( Figures S6B and S6C), with one exception. Anti-CTLA-4-naive and anti-CTLA-4-experienced tumors with BRAF V600 mutations that responded to anti-PD-1 had equally high TIARA-PD1 scores. The inverse relationship between TIDE and TIARA-PD1 was strongest in BRAF V600 and TWT tumors, and IMPRES did not stratify differences across clinical groups within any of the genetic subtypes ( Figure S6C). Together, the trends supported by TIDE and TIARA-PD-1 suggest that there may be differences in T cell infiltration and dysfunction across genetic subtypes, which could further implicate the concurrent or sequential targeting of these alterations by MAPK inhibitors in modulating the TME and response to anti-PD-1 therapy. 33 It is important to note that these patients were not annotated for their prior treatment with MAPK inhibitors.
Statistical models to predict response to anti-PD-1 in cutaneous melanoma are improved by including prior anti-CTLA-4 experience Differential gene expression analyses yielded consistent patterns across clinical subsets defined by prior anti-CTLA-4 experience and anti-PD-1 response; however, we aimed to identify which features of the defined genomic, transcriptomic, or clinical features would best predict anti-PD-1 response in cutaneous melanoma. Utilizing a LASSO-regularized logistic regression (LRLR) for variable selection, [34][35][36] we compared anti-PD-1responsive and nonresponsive tumors with WES and RNA-seq data from baseline tumors (n = 90). Features included statistically significant results across both data types from the previous analyses as well as published features (Table S3; Figure S7A). Features selected by the model that positively predicted anti-PD-1 response included PIK3C2G mutations, B cell signatures, and TIARA-PD-1, while features selected to predict anti-PD-1 nonresponse included expression of AREG or MAPK10 (Figures 6A and S7B). Some features were only selected based on the status of prior anti-CTLA-4 therapy, further demonstrating the interaction between prior anti-CTLA-4 therapy with gene expression features in predicting anti-PD-1 response, such as the Hallmark genes downregulated by KRAS activation as a predictor of anti-PD-1 nonresponse in anti-CTLA-4-naive tumors.
Because of the lack of genomic correlates in the integrated analysis and to increase statistical power, we attempted three additional approaches for all anti-PD-1-treated cutaneous melanoma tumors with RNA-seq data available (n = 131), regardless of matched WES. We performed LRLR with two distinct models. Model A was a complete model that integrates statistically significant factors related to response to anti-PD-1 regardless of association with anti-CTLA-4 experience, and model B only considered statistically significant factors without interaction with anti-CTLA-4 experience (Figures 6B, 6C, and S7C). Model A was subsequently divided into two approaches, utilizing either LRLR-selected features in a complete model containing interaction effects with anti-CTLA-4 experience (A1) or excluding those interaction effects (A2).
Features selected by these models, consisting of only transcriptomic features, were consistent with the first approach, which integrated genomic and transcriptomic features ( Figure 6A). However, there were additional features selected, including genes related to the extracellular matrix across all samples and the interaction effects in association with genes regulated by KRAS activation or MYC, myogenesis, and TIARA-PD-1 (Figures 6B and 6C). Exclusion of prior anti-CTLA-4 interactions in models A2 and B selected similar features as model A1, but inclusion of the interactions improved the leave-one-out cross-validation (LOOCV)-predicted probabilities and accuracy of the complete model ( Figures 6B  and 6C). Model A1 showed significant improvement in predicting anti-PD-1 response over models A2 (Wilcoxon rank-sum test, p = 8.1eÀ05) and B (Wilcoxon rank-sum test, p = 4.8eÀ07). This was shown by predicting the greatest probability of response in anti-PD-1 responding biopsies and the lowest probability of response in anti-PD-1 non-responding biopsies ( Figure 6C). Furthermore, model A1 had the highest overall accuracy in predicting anti-PD-1 response (74.8%), followed by A2 (65.6%) and B (62.6%) (Figure 6C). A random forest classifier selected features similar to our LRLR results with an accuracy of 66.7% (Figures S7D and  S7E). Collectively, this analysis suggests that annotation of anti-CTLA-4 experience is highly impactful in the context of studying anti-PD-1 response in cutaneous melanoma, and its inclusion may significantly refine the discovery of predictors or mechanisms responsible for anti-PD-1 response.

DISCUSSION
Correlates of ICI response have been studied using genomic, transcriptomic, and pathologic approaches in a global attempt to stratify patient populations through analysis of individual and aggregate clinical cohorts. 7,12 Gaining access to such datasets is a challenge in itself, and pooling larger cohorts requires additional methodology and expertise, scientific and clinic. Here, we harmonized the molecular and clinical annotation of seven clinical cohorts of patients with melanoma treated with ICIs, including two previously unpublished datasets. Integration of these studies allowed us to leverage the increased statistical power of an aggregate resource to more rigorously address the evaluation of anti-PD-1 response in cutaneous melanoma. We found that genomic and transcriptomic features predictive of anti-PD-1 response were modified by prior anti-CTLA-4 administration, highlighting key biological differences across baseline melanoma tumors studied across retrospective immunotherapy trials.
One of the primary benefits of harmonizing molecular datasets is the increased statistical power for more rigorous identification of biological features of patient subgroups. Other studies have attempted to increase statistical power by combining various clinical datasets, including those spanning different tumor types or treatment modalities. Alternatively, we focused on the utility of deeper clinical annotation, particularly with regards to melanoma subtype, treatment regimen, and treatment history. While Liu et al. 8 previously reported some correlates of anti-PD-1 response, controlling for melanoma subtype and anti-CTLA-4 experience, some were not reported as significant because of multiple testing correction. However, the rigorous harmonization of cohorts as well as addition of the CheckMate 064 and 067 datasets made it possible to detect signals associated with anti-PD-1 response, such as PIK3C2G missense mutations. This study highlights correlative findings that may support or drive experimental studies moving forward so that the mechanistic understanding of these correlations can be properly implemented clinically. Importantly, we found that refined stratification of cutaneous melanoma tumors by anti-CTLA-4 experience and anti-PD-1 response led to distinct sets of findings regarding

OPEN ACCESS
Article TMB, inflammatory gene expression patterns, and potential tumor-intrinsic gene expression patterns. The clinical datasets interrogated in this study comprised samples spanning melanoma subtypes; however, our assessment of TMB, confirming known differences in histologically defined subtypes, demonstrated the need for histopathologic annotation of samples for more rigorous analysis of molecular and clinical correlates. Collectively, this emphasizes why these efforts of data curation and harmonization are critical to clarify our understanding of predictive biomarkers for ICI treatment strategies and may also explain the lack of clarity, inconsistency, or lack of statistical signal in previous reports across these settings.
It is important to consider the relevance of our findings from these retrospective melanoma datasets in the context of the current treatment paradigm. Given the establishment of anti-PD-1 alone or in combination with anti-CTLA-4, as well as the recent approval of combined anti-PD-1/anti-LAG-3 therapy in the firstline setting for advanced melanoma, patients do not receive single-agent anti-CTLA-4 as a frontline therapy. Our analysis shows far fewer differences between anti-PD-1-responsive and nonresponsive tumors in the anti-CTLA-4-naive setting, which may highlight the need for refinement of predictive biomarkers to anti-PD-1, specifically in this setting. The known mechanism of PD-1 blockade involves recognition of tumor cells by CD8 T cells, and detection of infiltrating T cells in baseline samples has been a strong predictor of clinical benefit following anti-PD-1 therapy. 11 Anti-CTLA-4 stratified analysis of bulk RNA-seq tumor infiltration signatures revealed that prior CTLA-4 blockade was associated with significantly higher global immune infiltration, with enrichment for CD8 + T cell, B cell, TLS, and cDC1 signatures, which have all been implicated previously in anti-PD-1 responsiveness. 9,11 These results agree with pathology analysis data supporting that CTLA-4 blockade therapy increases intratumor T cell infiltrates regardless of clinical response to therapy. 29,37 There was also a lower enrichment of cell types commonly associated with immune suppression, such as neutrophil and Treg cell signatures, among patients with clinical benefit to anti-PD1. Furthermore, there were much lower signals associated with immune cell populations in anti-PD-1-nonresponsive tumors, regardless of anti-CTLA-4 experience.
We established TIARA-PD-1 as a summary metric of immune signatures differentially associated with response to anti-PD1 in cutaneous melanoma and found diverse transcriptional signatures as a function of immune infiltration based on anti-PD-1 response and prior anti-CTLA-4 experience. Specifically, an increase in TIARA-PD1 was strongly associated with a decrease in extracellular matrix, matrisome-associated proteins, glycoproteins, and myogenesis as well as glycolysis and hypoxia signatures. These signatures were consistent and more striking in anti-CTLA-4-experienced patients, while non-responders to anti-PD1 did not demonstrate a significant negative enrichment of many extracellular matrix (ECM) and hypoxia signatures despite an increase in TIARA-PD1. This is consistent with reports recognizing the ECM as a crucial microenvironmental component affecting the immune response in tumors as well as the role of a hypoxic TME in inducing immune suppression and resistance. [38][39][40][41] In addition, TIARA-PD-1 was more similar across anti-CTLA-4-experienced and -naive tumors specifically in BRAF V600-mutant tu-mors. These results warrant further study of combinatorial treatment strategies to enhance immunotherapy and potentiate cancer immunity for tumors resistant to immunotherapy despite strong immune infiltration. This should include immunotherapeutic strategies and targeted therapies, given the potential synergy in combining or sequencing MAPK inhibitor therapies in the context of ICIs. 42,43 Early ICI genomics studies reported the association between clinical response and TMB, 17,25,26,44 which has subsequently been used to positively identify ICI responders across tumor types. 7,12 Generally, it is thought that the presence of more mutations increases the chances of neoantigens arising that can be targeted by the immune system. However, in cancers that commonly have high TMB, such as melanoma, it can be difficult to identify specific recurrently mutated genes associated with response because there are many mutation events to consider when performing statistical tests. Our approach, assessing genes with corresponding mutation types in a larger dataset, enabled us to identify the positive association between PIK3C2G missense mutations and ICI response. However, because PIK3C2G is infrequently expressed overall in melanoma, the association between genomic and gene expression profiles remains unclear. 12 This is further supported by the lack of association between TMB and TIARA-PD-1, an aggregate metric of immune cell populations strongly associated with predicted anti-PD-1 response in cutaneous melanoma. Together, these findings suggest that defining the immunogenicity of tumors may require more specific or alternative genomic metrics in conjunction with assessment of immune recognition at the tumor site.
The reported, harmonized dataset presented in this study sets a foundation for investigating the heterogeneity of clinical responses to ICIs in melanoma. These harmonized molecular and clinical annotation efforts enabled us to focus on specific questions related to anti-PD-1 response in cutaneous melanoma, with the appropriate account of critical clinical and technical batch effects. These systematic approaches are necessary to advance the rigorous and complex analyses required to optimize ICI treatment strategies.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following: Tumor immunogenicity associated with response to anti-PD-1 (TIARA-PD-1) B EPIC, TIDE, and IMPRES B Model generation and evaluation d QUANTIFICATION AND STATISTICAL ANALYSIS Pieris, Bioentre, Iovance, Catalym, and Amgen; patents for methods for treating MHC class I polypeptide-related sequence A disorders (20100111973, pending, with royalties paid), tumor antigens and uses thereof (7250291, issued), angiopoiten-2 biomarkers predictive of anti-immune checkpoint response (20170248603, pending), compositions and methods for the identification, assessment, prevention, and treatment of melanoma using PD-L1 isoforms (20160340407, pending), therapeutic peptides (20160046716, 20140004112, 20170022275, and 20170008962, all pending, and 9402905, issued), methods of using pembrolizumab and trebananib (pending), vaccine compositions and methods for restoring NKG2D pathway function against cancers (10279021, issued), antibodies that bind to MHC class I polypeptide-related sequence A (10106611, issued), and anti-galectin antibody biomarkers predictive of anti-immune checkpoint and anti-angiogenesis responses (20170343552, pending); data safety monitoring board and advisory board participation for Aduro and Checkpoint Therapeutics; scientific advisory board leadership for Bicara and Apricity; and stock options in Checkpoint Therapeutics, Pionyr, Apricity, and Bicara. S.B., L.S., D.T., and T.T. are employees of BMS. C.N.S. received consulting fees from the Rare Cancer Research Foundation. D.K.W. is an employee, founder, and equity holder at Santa Ana Bio; holds equity in Immunai, and has received consulting fees from Rubius Therapeutics, DeepMind, Illumina, and Guardant Health. A.R. has received honoraria from consulting with Amgen, BMS, Chugai, Genentech, Merck, Novartis, Roche, and Sanofi; is or has been a member of the scientific advisory board and holds stock in Advaxis, Arcus Biosciences, Bioncotech Therapeutics, Compugen, CytomX, Five Prim, FLX-Bio, ImaginAb, Isoplexis, Kite-Gilead, Lutris Pharma, Merus, PACT Pharma, Rgenix, and Tango Therapeutics; and has received research funding from Agilent Technologies and BMS through Stand Up to Cancer.