Distinct Mechanisms of Mismatch-Repair Deficiency Delineate Two Modes of Response to Anti–PD-1 Immunotherapy in Endometrial Carcinoma

Mismatch repair–deficient (MMRd) cancers are characterized by exceptionally high mutational loads, particularly at repetitive microsatellite regions that are prone to replication errors [microsatellite instability-high (MSI-H); refs. 1, 2]. The high mutation burden in MSI-H/MMRd tumors is thought to result in the generation of ample tumor-specific neoantigens that can subsequently be targeted by the immune system in the setting of immune-checkpoint blockade (ICB; refs. 3, 4). As the first tissue-agnostic ICB therapy approved by the FDA, the PD-1 checkpoint inhibitor pembrolizumab received accelerated approval in 2017 for the treatment of MSI-H/MMRd solid tumors that have progressed on prior lines of therapy (5–8). However, even among MSI-H/MMRd tumors, there is considerable variability in response rates to ICB, both within and across cancer types. For patients with MSI-H/MMRd colorectal cancer or endometrial cancer, 44% and 57% of patients, respectively, demonstrated an objective radiographic response to anti–PD-1 ICB (5). In stark contrast, only 18% of patients with MSI-H/MMRd pancreatic cancer and 0% of patients with brain cancer responded to anti–PD-1 ICB (7). To increase the proportion of patients who derive benefit from ICB, it is crucial to identify the genetic, immunologic, and environmental factors underlying these variations in response.

MSI-H/MMRd status is commonly determined by PCR of microsatellite regions and/or IHC for key MMR proteins (9), and such assays were utilized in the clinical trials that demonstrated the efficacy of pembrolizumab for MSI-H/MMRd solid tumors. However, neither PCR nor IHC can provide insight into the underlying molecular mechanisms of MMRd. Whereas some tumors exhibit somatic loss-of-function mutations in key MMR genes (hereafter termed mutational MMRd, or mut-MMRd), more than 70% of MSI-H/MMRd endometrial cancer cases (10) instead exhibit epigenetic silencing of the MLH1 gene promoter through DNA methylation (epigenetic MMRd, or epi-MMRd; ref. 9).

We recently conducted a phase II trial (NCT02899793) evaluating the efficacy of pembrolizumab in patients with MSI-H/MMRd endometrial cancers (n = 24 evaluable patients; ref. 11). In this trial, we classified patients into two groups based on the putative molecular mechanism driving MMRd (i.e., mut-MMRd or epi-MMRd). We observed a 100% objective response rate (ORR) for mut-MMRd patients, compared with an ORR of 44% for epi-MMRd patients (11). Accordingly, mut-MMRd patients had a more rapid clinical response, as well as significantly longer progression-free survival and overall survival compared with epi-MMRd patients, suggesting important differences in the immunologic response to anti–PD-1 ICB across these two subtypes of MSI-H/MMRd tumors.

We therefore hypothesized that certain immunologic parameters may vary between patients with mut-MMRd or epi-MMRd tumors, providing a potential mechanistic explanation for the variable therapeutic benefit of anti–PD-1 ICB across MSI-H/MMRd endometrial cancers. Here, we report a comprehensive characterization of the tumor genomes and circulating immune responses in the participants of our recent phase II clinical trial. We performed single-cell RNA sequencing (scRNA-seq) and matched T-cell receptor repertoire sequencing (scTCR-seq) of peripheral blood mononuclear cells (PBMC) from these patients, before and after treatment with pembrolizumab, along with pretreatment bulk TCR-seq of tumor-infiltrating lymphocytes (TIL) and tumor exome sequencing. By interrogating this data set, we uncovered two distinct axes of antitumor immunity that appear to be differentially engaged by anti–PD-1 ICB in patients with mut-MMRd versus epi-MMRd tumors.

Study participants were patients with recurrent MSI-H/MMRd endometrial carcinoma who enrolled in our previously described single-arm, open-label phase II clinical trial evaluating the efficacy of pembrolizumab in this population (n = 24 evaluable patients; ref. 11). The molecular mechanism of MMRd was determined by MLH1 promoter methylation analysis, FoundationOne profiling, whole-exome sequencing (WES), PCR, and IHC for MMR proteins. Six patients were classified as mut-MMRd and the remaining 18 patients as epi-MMRd (Supplementary Table S1; ref. 11).

Response to pembrolizumab (200 mg of pembrolizumab every 3 weeks for ≤2 years) was assessed by RECIST 1.1 criteria (12), defining responders as patients with a ≥30% decrease in tumor size (Fig. 1A). All mut-MMRd patients demonstrated an objective response to pembrolizumab (6/6, 100%), whereas 44% of epi-MMRd patients responded (8/18; Fig. 1B and C). As we were interested in understanding the factors that distinguished epi-MMRd responders from nonresponders, we subsequently categorized the patients into 3 groups for further analysis: epi-MMRd nonresponders (NR; n = 10), epi-MMRd responders (epiR; n = 8), and mut-MMRd responders (mutR; n = 6). Patient ages were similar across the three groups (Supplementary Fig. S1A).

Longitudinal analysis of tumor sizes revealed that responses to pembrolizumab were sustained for well over a year, with faster response kinetics in mutR patients (Fig. 1D). As we previously reported, mut-MMRd patients had longer progression-free survival and overall survival compared with epi-MMRd patients (Fig. 1E; Supplementary Fig. S1B and S1C). Tumor stage and grade were comparable across the three groups of patients (Fig. 1F and G). Importantly, tumor mutation burden (TMB) was significantly higher in mutR patients compared with NR or epiR patients, whereas TMB was not significantly different between NR and epiR patients (Fig. 1H). Using a set of high-confidence somatic exonic mutations, we predicted the number of MHC-I and MHC-II neoantigens in each sample, revealing a similar trend toward increased neoantigen counts in mutR patients, though not statistically significant (Supplementary Fig. S2A–S2C; Supplementary Table S2).

To corroborate these findings in an independent cohort, we examined the uterine endometrial carcinoma cohort from The Cancer Genome Atlas (TCGA UCEC; ref. 13). In this cohort, 75.4% of the tumors were classified as endometrioid carcinomas. We categorized the tumors into four groups: (i) MMR-proficient (MMRp), (ii) epi-MMRd, (iii) mut-MMRd, and (iv) MMRd, not otherwise specified. To define epi-MMRd tumors, we identified tumors that were previously classified as MSI-H and had a methylation β-value ≥0.5 for the MLH1 promoter (Supplementary Fig. S2D). To define mut-MMRd samples, we identified tumors that were previously classified as MSI-H (14) and that exhibited mutations in a canonical MMR gene (MLH1, MSH2, MSH6, or PMS2), but did not have MLH1 promoter methylation. As a validation of the β-value threshold, we confirmed that samples with a methylation beta-value ≥0.5 for the MLH1 promoter had significantly decreased expression of MLH1 (Supplementary Fig. S2E). There was no association between tumor stage and the four categories of tumors (Supplementary Fig. S2F). Comparing across groups, we observed that all MSI-H/MMRd groups had significantly higher TMB than MMRp tumors, as expected (Fig. 1I). We further observed that mut-MMRd tumors had significantly higher mutation rates than epi-MMRd tumors, confirming our observations with the trial cohort (Fig. 1I). Similarly, all MSI-H/MMRd groups had significantly higher predicted neoantigen load (14) than MMRp tumors (Fig. 1J), with mut-MMRd tumors exhibiting the highest numbers of predicted neoantigens. Thus, analysis of an independent endometrial cancer cohort confirmed the relationship between the molecular mechanism of MMRd and TMB.

Of note, the patients profiled in the TCGA UCEC cohort were not treated with ICB therapy. To further investigate whether mut-MMRd tumors have higher ICB response rates than epi-MMRd tumors, we reanalyzed data from the MSK-IMPACT cohort of cancer patients treated with ICB therapy (15–17). Because methylation data were not available for this cohort, we instead categorized the patients with endometrial cancer into three groups: (i) MMRp, (ii) mut-MMRd, and (iii) nonmutational MMRd (nonmut-MMRd). In the MSK-IMPACT cohort, mut-MMRd patients were defined by the combination of MSI-H status and the presence of mutations in MLH1, MSH2, MSH6, or PMS2. In contrast, nonmut-MMRd patients were defined by the combination of MSI-H status and the absence of mutations in MLH1, MSH2, MSH6, or PMS2. The classification scheme for nonmut-MMRd tumors in the MSK-IMPACT cohort was thus different from that used to define epi-MMRd in our own cohort or in the TCGA cohort; this approach represented the closest possible surrogate given the available data. Comparing the response rates with ICB therapy across the three groups, we observed that mut-MMRd patients had the highest response rate (5/8, or 63%), whereas nonmut-MMRd patients and MMRp patients had response rates of 36% and 23%, respectively (Fig. 1K). These data further corroborate our finding that mut-MMRd tumors have higher response rates to ICB therapy than nonmutational MMRd tumors, albeit with the limitation that MLH1 promoter methylation data were not available for the MSK-IMPACT cohort.

Given the association between mut-MMRd and higher TMB, one might expect that differences in tumor immunogenicity due to neoantigen load are sufficient to explain the varying ICB responses between mut-MMRd and epi-MMRd patients (4). However, when we further stratified patients by response status, we observed no significant differences in TMB between responders and nonresponders within each group (Fig. 1L). This finding confirmed our earlier observation that NR and epiR patients had comparable TMB (Fig. 1H). Thus, although the molecular mechanism of MMRd is reproducibly associated with TMB, in the setting of elevated neoantigen load (e.g., among MSI-H/MMRd endometrial cancers), factors beyond TMB are crucial in shaping the antitumor immune response upon ICB therapy.

Analysis of computationally predicted pathogenic and deleterious mutations revealed several recurrently mutated genes across the patient cohort (Fig. 2A), including PTEN, ARID1A, RPL22, KMT2B, and KMT2D (Supplementary Table S3). Across the entire cohort, only mutations in ARID1A and CTCF were nominally associated with response to pembrolizumab, with the limitation of the small sample size for this clinical study of phase II trial participants (Fig. 2B). CTCF encodes a key factor for three-dimensional (3D) chromatin organization (18), and CTCF binding sites are often mutated in cancer (19). We observed that mutations in CTCF were associated with response to pembrolizumab (Fig. 2C), suggesting that dysregulated 3D chromatin structure may promote immunotherapy response in this context. On the other hand, ARID1A encodes a subunit of the BAF SWI/SNF chromatin remodeling complex (20) and is frequently mutated across multiple cancer types (21). Although ARID1A deficiency has previously been associated with enhanced ICB responsiveness (22), possibly through loss of its interactions with MMR machinery, we found that ARID1A mutations were associated with resistance to pembrolizumab (Fig. 2C). This difference may stem from the fact that our current study was restricted to patients with MSI-H/MMRd tumors, as the sensitizing effect of ARID1A deficiency through loss of MMR function would not be relevant in this setting.

Prior studies across multiple contexts have indicated that JAK1 mutations can protect tumor cells from CD8⁺ T-cell cytotoxicity by disabling IFNγ pathway signal transduction (23–31). For instance, multiple CRISPR screens using in vitro coculture systems have identified JAK1 loss as a key mechanism by which tumor cells can evade T-cell cytotoxicity (25, 32–34). Analogous observations have been described in the clinical setting. For instance, a case series of melanomas with acquired resistance to pembrolizumab demonstrated that new JAK1, JAK2, or B2M mutations can arise following an initial response to ICB (24). Similarly, other case series have described the potential impact of pretreatment JAK1 mutations on initial ICB responses, concluding that JAK1 loss can also confer primary resistance to ICB (23, 35). Although these reports are compelling, the relatively small sample sizes of these prior clinical studies are an important limitation. To ascertain whether pretreatment JAK1 mutations indeed confer resistance to ICB, larger cohorts of JAK1-mutant and JAK1-nonmutant tumors are needed.

Our cohort of patients with MSI-H/MMRd endometrial cancer provided an opportunity to explore this question. Of the tumors profiled in this study, 42% (10 of 24) exhibited predicted pathogenic/deleterious mutations in JAK1 prior to treatment (Fig. 2A). Analysis of the TCGA Pan-Cancer data set (36) revealed the highest rates of JAK1 alterations (14%) in the UCEC cohort (Supplementary Fig. S3A), further increasing to 30% of MSI-H endometrial cancers (37). Analysis of the AACR-GENIE cohort (38) showed similar findings, with the highest rates of JAK1 mutations (11%) in endometrial cancers (Supplementary Fig. S3B). These data indicate that JAK1 mutations are abundant in endometrial cancers, possibly hinting at active immune surveillance even in the absence of ICB therapy.

However, contrary to our expectations, pretreatment JAK1 mutations were not associated with primary resistance to pembrolizumab in this cohort of MSI-H/MMRd endometrial cancers (Fig. 2B). In fact, 70% (7 of 10) of the JAK1-mutant tumors responded to pembrolizumab, compared with 50% (7 of 14) of the nonmutant tumors (Fig. 2C). To understand the clonality of these JAK1 mutations, we converted the variant allele frequencies (VAF) into cancer cell fractions (CCF). By accounting for tumor cellularity and allele-specific copy-number variation, CCFs more accurately represent the clonality of mutations compared with VAFs (ref. 39; Methods). We observed that JAK1-mutant CCFs were consistently above 0.25 in JAK1-mutant tumors, indicating that JAK1 mutations were present in more than 25% of the cells in each tumor (Fig. 2D; Supplementary Fig. S3C and S3D). Of note, all three of the patients who exhibited a complete response to pembrolizumab (PEM06, PEM23, and PEM25) had JAK1 mutations prior to treatment (Fig. 2A). In these patients, JAK1-mutant cells were estimated to comprise a substantial proportion (>80%) of the total tumor cell population (Fig. 2D). If JAK1 mutations indeed conferred primary resistance to immunotherapy, these JAK1-mutant cells would be expected to drive disease progression while on pembrolizumab. Thus, it is notable that a complete response was observed in patients PEM06, PEM23, and PEM25, even though their tumors were largely comprised of JAK1-mutant cells.

To further understand the relationship between JAK1 alterations and ICB response, we performed WES on 8 epi-MMRd tumors after the initiation of pembrolizumab. We observed that three of the tumors (two NR and one epiR) had JAK1 frameshift mutations that were not detected in the pretreatment samples (Fig. 2D; Supplementary Fig. S3D). These JAK1 mutations could represent clonal selection in the setting of enhanced antitumor immune pressure following pembrolizumab, consistent with acquired resistance to ICB. We speculate that rare JAK1-mutant subclones were already present prior to treatment initiation, but given the limitations of single-site exome sequencing for capturing infrequent subclonal variants in heterogeneous tumors, these JAK1 variants were not detected pretreatment.

On the other hand, one epiR patient (PEM19) had a JAK1 frameshift mutation (p.M316fs; c.947_948insAC) at baseline that subsequently was not detected on disease progression while on pembrolizumab (Supplementary Fig. S3D). This suggests that the pretreatment JAK1 mutation in patient PEM19 was not under positive selection despite the heightened immune pressure from ICB treatment. Furthermore, it indicates that the pretreatment JAK1-mutant subclone was unlikely to be responsible for driving disease progression, as an ICB-resistant subclone would be expected to constitute a large enough proportion of the posttreatment tumor to be detectable on WES. With that being said, a caveat to this interpretation is that the JAK1 mutation may not have been identified in the posttreatment sample due to tumor heterogeneity and the limitations of single-site exome sequencing. Another possible explanation for the loss of this JAK1-mutant subclone is that an immune response was mounted against a neoantigen encoded by the JAK1 mutation. The top predicted mutant MHC-I epitope at this site showed low predicted affinity (IC₅₀: 5,434.3 nmol/L), whereas the top predicted MHC-II epitope had a moderate predicted affinity (354.46 nmol/L). It is therefore possible that tumor cells carrying this JAK1 mutation were eliminated in patient PEM19 following ICB therapy due to a T-cell response against this neoepitope.

Taken together, these data indicate that JAK1 mutations may be advantageous for endometrial cancer outgrowth by subverting immune surveillance. This was especially apparent in the case of three patients whose tumors appeared to have acquired JAK1 mutations following pembrolizumab initiation. Interestingly, however, we also found that pretreatment JAK1 mutations did not necessarily confer primary resistance to pembrolizumab. For instance, all three tumors with a complete response to anti–PD-1 ICB had near-clonal JAK1 mutations at baseline.

We were intrigued that only a fraction of patients with epi-MMRd tumors showed a response to anti–PD-1 immunotherapy, whereas all patients with mut-MMRd tumors had responded (Fig. 1B and C). Although mut-MMRd tumors had higher TMB than epi-MMRd tumors, TMB alone could not explain the observed variation in ICB response (Fig. 1H and L). We thus sought to explore the immunologic mechanisms driving differential ICB responses across patients. To do so, we sampled PBMCs before and after pembrolizumab treatment (n = 22 patients). In total, we performed paired scRNA-seq and scTCR-seq on more than 260,000 PBMCs (n = 52 samples; Supplementary Fig. S4A). We also performed bulk TCR-seq on the pretreatment TILs to construct a reference library of tumor-reactive TCR sequences.

After standardized data preprocessing, we visualized the scRNA-seq data through uniform manifold approximation and projection (UMAP; Fig. 3A). The cells were classified into cell types by analysis of differentially expressed genes (DEG) within each cluster (Supplementary Fig. S5A). We then further subclustered the T-cell populations for a more detailed perspective on the functional status of the circulating T cells in these patients (Fig. 3B; Supplementary Fig. S5B). We ultimately classified the PBMCs into 21 cell types (Supplementary Table S4). After filtering out samples with fewer than 1,000 total cells (n = 3 samples removed; Supplementary Fig. S4A), we compared the relative frequencies of the 21 cell types across time points and patient categories. At the level of the 21 cell clusters, we did not observe major systematic differences in cell type frequencies between NR, epiR, and mutR patients (Supplementary Fig. S5C). Similarly, there were inconsistent changes in cell type frequencies before and after pembrolizumab treatment, comparing within each category of patients (Supplementary Fig. S5C).

We reasoned that finer-resolution changes in specific cell subsets may be masked when analyzing the relative frequencies of broad cluster-based cell-type classifications. To identify subtler shifts in circulating immune cell subpopulations and cell states, we grouped the cells into “neighborhoods” through k-nearest neighbor clustering, as implemented through the Milo analysis framework (40). Differential abundance analysis of cell neighborhoods revealed that prior to treatment, epiR patients did not have significant enrichment for any specific subpopulations compared with NR patients (Supplementary Fig. S6A). In contrast, mutR patients exhibited significant enrichment for specific subpopulations of activated CD8⁺ T cells when compared with NR patients (Fig. 3C). DEG analysis showed that the specific CD8⁺ T-cell neighborhoods enriched in mutR patients were defined by high expression of KLRG1, EOMES, and LAG3, along with lower expression of inhibitory receptors PDCD1, CTLA4, and KLRC1 (encoding NKG2A; ref. 41; Fig. 3D). There was no difference in IFNG or TNF expression between enriched and nonenriched CD8⁺ T-cell neighborhoods (Supplementary Fig. S6B). In aggregate, this expression signature suggests that a CD8⁺ T-cell effector population is poised for mounting antitumor responses in mutR patients.

After treatment with anti–PD-1 ICB, epiR patients again did not show any significant differences in cell neighborhood abundances (Supplementary Fig. S6C), whereas mutR patients were enriched for several CD8⁺ T-cell neighborhoods (Fig. 3E). The mutR-enriched CD8⁺ T-cell neighborhoods in the posttreatment setting were characterized by elevated expression of several cytolytic effector molecules (GZMB, GNLY, and TNF; Fig. 3F) but not IFNG (Supplementary Fig. S6D). CD8⁺ T neighborhoods enriched posttreatment in mutR patients also demonstrated reduced expression of inhibitory receptors such as HAVCR2 (encoding TIM3), PDCD1, and LAG3. Overall, these data point to the presence of an activated CD8⁺ T-cell population in mutR patients that is associated with antitumor responses following anti–PD-1 ICB.

To further characterize the T-cell responses in these patients, we next focused on the scTCR-seq repertoires. Although there were no statistically significant differences in TCR clonality (by the Simpson clonality index) among the three categories of patients before or after treatment, there was a trend toward increased TCR clonality in responders, particularly in a subset of mutR patients (Supplementary Fig. S6E and S6F). Using the pretreatment bulk TCR-seq data collected from TILs, we generated a reference library of tumor-reactive TCR sequences in each patient, which we used to identify putative tumor-reactive clonotypes among the peripheral T cells. With this reference, we quantified the percentage of circulating T cells with TIL-matching TCR sequences in each sample. Prior to treatment, there were no significant differences in the relative proportions of TIL-matched circulating T cells among the different patient groups (Supplementary Fig. S6G). In contrast, there was a significant increase in the percentage of TIL-matched circulating T cells in mutR patients compared with NR (Supplementary Fig. S6H). Of note, the expansion of peripheral T-cell populations has previously been associated with tumor infiltration and immunotherapy response across cancer types (42). Collectively, these data revealed a circulating effector CD8⁺ T-cell response in mutR patients that appeared to be further amplified by anti–PD-1 immunotherapy.

Our earlier analysis of the tumor exome sequencing data revealed that pretreatment JAK1 alterations were not associated with primary resistance to pembrolizumab (Fig. 2B and C). We therefore wondered if patients with JAK1-mutant tumors had different circulating immune responses compared with those with JAK1-nonmutant tumors. We performed Milo analyses comparing JAK1-mutant and JAK1- nonmutant tumors within each patient group (NR, epiR, and mutR). This analysis revealed minimal changes in cell neighborhood abundances when compared by JAK1 mutation status (Supplementary Fig. S7A and S7B), suggesting that patients with JAK1-mutant tumors did not have systematic differences in the composition of the circulating immune response. However, we cannot entirely rule out that immune responses against JAK1-mutant tumors may subtly differ from those against JAK1-nonmutant tumors, as the present analysis is limited by the relatively small size of the trial cohort.

Distinct from mutR patients, epiR patients had similar TMB to NR patients (Fig. 1H) and did not exhibit significant enrichment for activated CD8⁺ T cells (Supplementary Fig. S6A and S6B). These data suggest that the underlying mechanism of MMRd (i.e., mutational or epigenetic) influences the intensity of MSI, potentially leading to differences in the intrinsic characteristics of the associated tumors and the cellular composition of the circulating immune response. Although these factors could explain the high response rates against mut-MMRd tumors, they could not explain why approximately half of the patients with epi-MMRd tumors responded to anti–PD-1 ICB. We therefore hypothesized that transcriptional characteristics of the circulating immune response could distinguish patients with epi-MMRd tumors that respond to PD-1 blockade from those that do not.

We subsequently set out to identify the transcriptional features of T and natural killer (NK) cells in epiR patients compared with NR patients. However, as the number of profiled cells varied widely between samples (Supplementary Fig. S4A), conventional scRNA-seq approaches for differential expression such as the Mann–Whitney test (the default in the widely used Seurat analysis package; refs. 43, 44) would be weighted toward samples with more cells. In a similar vein, we were concerned that considering each single cell as an independent event would artificially inflate statistical power. We instead used a random sampling–based approach to partition the cells in each sample into “pseudobulk” profiles (see Methods). Specifically, we generated five pseudobulk profiles for each cell type in each sample, then performed DEG analysis using a generalized linear model as implemented in DESeq2 (45). This approach equally weights each sample and avoids overcounting single cells as independent events, providing advantages over the Mann–Whitney test when analyzing large multisubject scRNA-seq data sets.

Comparing responders to nonresponders prior to anti–PD-1 immunotherapy, we quantified the number of significant DEGs across T and NK cell populations (Supplementary Fig. S8A). Pathway analysis revealed that CD16⁺ NK cells in epiR patients demonstrated strong upregulation of genes involved in the proteasome, NK cell–mediated cytotoxicity, oxidative phosphorylation, and NOD-like receptor signaling (Supplementary Fig. S8B), whereas these changes were not seen in CD16⁺ NK cells from mutR patients. Across these T or NK cell populations, the majority of upregulated pathways were identified in CD16⁺ NK cells. These data indicate that prior to anti–PD-1 immunotherapy initiation, CD16⁺ NK cells in epiR patients were transcriptionally distinct from those in NR patients, with higher expression of genes that are relevant to antitumor immunity.

We observed a similar pattern in the posttreatment setting, comparing epiR patients and mutR patients with nonresponders (Supplementary Fig. S8C). After pembrolizumab treatment, activated CD8⁺ T cells in epiR patients uniquely demonstrated upregulation of genes involved in oxidative phosphorylation and the proteasome (Supplementary Fig. S8D). In a similar manner, activated CD4⁺ T cells showed enrichment for NOD-like receptor signaling and fatty acid metabolism in epiR patients but not in mutR patients. On the other hand, the predominant significant pathway-level changes in NK cells after treatment were mutR-specific downregulation events. For instance, CD16⁺ NK cells and CD56⁺ NK cells in mutR patients both had significant downregulation of NK cell–mediated cytotoxicity, but these pathway-level alterations were not observed in epiR patients.

We next identified DEGs in a longitudinal manner within each patient group, comparing pretreatment and posttreatment transcriptional profiles for each cell type (Supplementary Fig. S9A). For instance, when comparing pre- and posttreatment samples, activated CD4⁺ T cells, CD8⁺ T cells, and CD16⁺ NK cells showed more upregulated and downregulated DEGs in epiR patients than NR or mutR patients. On the pathway level, both naïve and activated CD4⁺ T cells in epiR patients exhibited upregulation of antigen processing and presentation, along with downregulation of oxidative phosphorylation (Supplementary Fig. S9B). CD16⁺ NK cells in epiR patients uniquely demonstrated upregulation of genes involved in the phagosome, tight junctions, focal adhesion, and leukocyte transendothelial migration.

These analyses suggested differing modes of antitumor immunity in epiR patients compared with mutR patients. We thus directly compared the transcriptional profiles of immune cells from epiR and mutR patients, stratifying by pretreatment and posttreatment time points (Fig. 4A). Pathway analyses of the resulting DEGs showed relative enrichment of TCR signaling and TNF signaling in activated CD8⁺ T cells from mutR patients compared with epiR patients (Fig. 4B), corroborating our earlier observation that mutR patients are characterized by potent effector CD8⁺ T cells (Fig. 3C–F). Conversely, CD16⁺ NK cells from epiR patients had higher expression of genes involved in NK cell–mediated cytotoxicity, graft versus host disease, and NOD-like receptor signaling (Fig. 4B). The enrichment of these pathways in CD16⁺ NK cells suggests that NK cells may play a key role in mounting antitumor responses in epiR patients.

Collectively, these findings indicate that the transcriptional profiles of T and NK cell populations in epiR patients are distinct from those of NR and mutR patients, both before and after pembrolizumab treatment. Although epiR patients did not show significant numerical enrichment of CD8⁺ T cells over their nonresponding epi-MMRd counterparts, transcriptional profiling pointed to heightened functionality of circulating NK cells in epiR patients.

The DEG analyses pointed to enhanced functionality of NK cells in epiR patients compared with NR patients. Given that mutR patients, but not epiR patients, had increased TMB and activated CD8⁺ T-cell frequencies compared with NR patients (Fig. 3C–F), and that T cells from mutR patients had comparatively higher expression of TCR signaling–related genes (Fig. 4B), we wondered if NK cells might be responsible for promoting antitumor responses in epiR patients. This notion was further supported by the observation that CD16⁺ NK cells in epiR patients showed higher expression of genes involved in NK cell–mediated cytotoxicity compared with their mutR counterparts (Fig. 4B).

To explore this possibility, we first examined the TCGA UCEC cohort of patients with endometrial cancer to identify tumor features that correlate with overall survival (13). Although the TCGA UCEC study predates the advent of ICB therapies for this patient population, we reasoned that by understanding the features that associate with survival at baseline (i.e., in the absence of ICB therapy), we could nevertheless identify factors that are clinically relevant for endometrial cancer biology. We constructed a multivariable Cox proportional hazards model examining a number of factors that have previously been associated with response to immunotherapy (14, 46): nonsilent mutation rate (47), homologous recombination defects (48, 49), TGFβ response (50, 51), IFNγ response (52, 53), CD8⁺ T-cell abundance (47, 53–55), indel burden (56), MSI score (ref. 4; determined by MSIsensor, ref. 57), and NK cell abundance (58, 59). In the fully adjusted model, we found that homologous recombination defects were significantly associated with worse survival (HR, 1.39; 95% CI, 1.13–1.70), whereas activated NK cells were associated with better survival (HR, 0.72; 95% CI, 0.55–0.95; Fig. 4C). In contrast, CD8⁺ T cells were not significantly associated with survival in the multivariable analysis, though there was a trend toward improved survival (HR, 0.87; 95% CI, 0.66–1.15). Similarly, the elevated nonsilent mutation rate was associated with improved survival but did not reach statistical significance (HR, 0.64; 95% CI, 0.37–1.10).

The significant association between activated NK cells and overall survival suggests that NK cells might play an important role in immune responses against endometrial cancers. We therefore further interrogated the transcriptional characteristics of NK cells in epiR patients, compared with NR and mutR patients. We specifically focused on the CD16⁺ NK cell subset because CD16⁺ NK cells are the predominant NK cell subset in the circulation, have been demonstrated to be more naturally cytotoxic than CD56⁺ NK cells, and are preferentially recruited from the circulation to sites of inflammation (60, 61). In comparing the expression profiles of CD16⁺ NK cells from epiR patients to those of NR and mutR patients, we identified 248 genes that were commonly upregulated across both comparisons prior to ICB treatment (Fig. 4D). These included a number of genes relevant for NK cell functionality, such as ARPC3, CX3CR1, GZMA, LTB, and TNFSF10 (encoding TRAIL; refs. 62, 63). Gene ontology analysis showed significant enrichment for categories such as the exosome, ruffle membrane, TNF-mediated signaling, focal adhesion, and Rho cell motility signaling (Fig. 4E), indicative of an activated NK cell state.

Similarly, in the posttreatment setting, CD16⁺ NK cells from epiR patients had 290 commonly upregulated genes in comparison with NR and mutR patients (Supplementary Fig. S10A). These included genes that were already upregulated prior to treatment, such as ARPC3, CX3CR1, GZMA, and LTB, as well as posttreatment-specific DEGs, including EOMES, FLT3LG, LTA, and STAT1. Gene ontology analysis similarly revealed enrichment for the exosome and TNF-mediated signaling, along with upregulation of the proteasome, mitochondrion, immune response, and Fc-epsilon receptor signaling (Supplementary Fig. S10B).

As for downregulated DEGs, 131 genes (prior to treatment) and 363 genes (after treatment) were commonly downregulated in CD16⁺ NK cells from epiR patients compared with NR and mutR patients (Fig. 4F). Gene ontology analysis of the downregulated DEGs showed enrichment for ribonucleoprotein, nonsense-mediated decay, and poly(A) RNA-binding, both before and after anti–PD-1 ICB (Fig. 4G). We noted that ZFP36L2 was among the most strongly downregulated genes in CD16⁺ NK cells from epiR patients (Supplementary Fig. S10C). The ZFP36 family of RNA-binding proteins has previously been shown to translationally repress cytokine production in T cells by binding to 3′-untranslated regions, and loss of ZFP36 family members leads to enhanced antiviral responses (64, 65).

Given that activated NK cell abundances were associated with survival in the TCGA UCEC cohort (Fig. 4C), we wondered if the gene signatures of CD16⁺ NK cells in epiR patients could stratify patients with endometrial cancer into prognostic subgroups. We constructed a Lasso Cox regression model that included all 111 genes that were concordantly upregulated or downregulated in CD16⁺ NK cells from epiR patients, taking the intersection of pretreatment and posttreatment DEGs. We arrived at a final set of 4 genes (termed epiR-NK4) that were each significantly associated with survival in a multivariable analysis and that demonstrated a reciprocal relationship between survival and differential expression in epiR CD16⁺ NK cells (Fig. 4H). The epiR-NK4 gene set consisted of CD63, PPIB, CEBPB, and LDOC1, which were consistently upregulated (CD63 and PPIB) or downregulated (CEBPB and LDOC1) in CD16⁺ NK cells from epiR patients, both before and after PD-1 treatment. High expression levels of the epiR-upregulated genes CD63 and PPIB were associated with improved survival, whereas high expression levels of the epiR-downregulated genes CEBPB and LDOC1 were associated with worse survival (Fig. 4H).

To derive a summative statistic for the epiR-NK4 gene set across the TCGA UCEC cohort, we computed the first principal component (PC1) of the RNA-seq expression values for CD63, PPIB, CEBPB, and LDOC1 (Fig. 4I; Supplementary Fig. S10D). The PC1 variable loadings were positive for CD63 and PPIB, but negative for CEBPB and LDOC1; thus, samples with high expression of CD63/PPIB and low expression of CEBPB/LDOC1 would have a high epiR-NK4 score. We evaluated whether the epiR-NK4 score could stratify survival in specific subgroups of interest: (i) tumors with a high abundance of activated NK cells and (ii) tumors with MSI-H/POLE-hypermutated status. In both subgroup analyses, we found that the median epi-NK4 score defined two distinct prognostic categories, with a high epiR-NK4 score being associated with better survival (Fig. 4J and K). These findings suggest that the transcriptional signatures of CD16⁺ NK cells in epiR patients are associated with improved survival in patients with endometrial cancer.

A substantial proportion of MSI-H/MMRd tumors do not respond to therapy despite their exceptionally high mutational loads (7). We show here that among MSI-H/MMRd endometrial cancers, the molecular mechanism of MMRd is associated with response to anti–PD-1 ICB. As a potential explanation for why tumors with mut-MMRd have more robust responses compared with those with epi-MMRd, we found that the nature of the circulating immune response differs depending on the mechanism of MMRd. For patients with mut-MMRd tumors, the circulating immune repertoire is enriched for a subpopulation of effector CD8⁺ T cells; for patients with epi-MMRd tumors that respond to anti–PD-1 ICB, the circulating immune response is defined by a highly active CD16⁺ NK cell population.

We speculate that mut-MMRd tumors and epi-MMRd tumors have distinct tumor-intrinsic features that shape the nature of the circulating immune system and subsequent responses to anti–PD-1 ICB. Although mut-MMRd tumors exhibited higher TMB than epi-MMRd tumors, we also found that TMB could not fully explain the distinction between responders and nonresponders (Fig. 1H and L). Of note, prior studies have demonstrated that out of all the nonsynonymous mutations present within a tumor, only a small fraction will be recognized by TILs (66, 67). Considering that even a single neoantigen can be sufficient to drive potent antitumor immunity (68, 69), it is likely that the reported association between TMB and ICB response (16, 47) exhibits diminishing returns. Within a cohort of MSI-H/MMRd tumors that all have high neoantigen loads (as in the present study), neoantigen availability is less likely to be the key limiting factor for mounting effective antitumor immune responses. Thus, we hypothesize that other tumor-intrinsic differences beyond TMB collectively influence the state of the host immune system, leading to differential responses to anti–PD-1 ICB in mut-MMRd versus epi-MMRd tumors.

To that end, how might the molecular mechanism of MMRd influence the nature of the circulating antitumor immune response? Although answering this question will require further studies in controlled experimental systems, we hypothesize that there are at least two contributing factors.

First, the observed differences between mut-MMRd and epi-MMRd tumors likely stem in part from the rewritable nature of epigenetic alterations. Whereas all downstream progeny of a mut-MMRd tumor clone will carry the same MMR gene mutations, the extent of MMR silencing through MLH1 promoter methylation may vary within an epi-MMRd tumor (70). This heterogeneity in MLH1 silencing can occur at two separate levels: across different cells within a tumor and across time, for instance, as the microenvironmental state of the tumor evolves (71). To the extent that transient restoration of MLH1 expression (and thus MMR capacity) occurs in epi-MMRd tumors, it is plausible that a fraction of the mutations that had previously accumulated as a result of MLH1 promoter methylation will then be resolved (72, 73). An analogous phenomenon was recently described for BRCA1, where BRCA1-methylated tumors could adaptively regain BRCA1 function and acquire resistance to platinum therapies (74). This hypothesis is further supported by the observation that mut-MMRd tumors had significantly higher TMB than epi-MMRd tumors. As a consequence, some of the neoantigens that were being targeted by CD8⁺ T cells may be lost in a subset of epi-MMRd tumor cells. In this scenario, NK cells may play a larger role in antitumor responses, as they can efficiently mount cytotoxic responses in a neoantigen-independent manner.

Second, epi-MMRd tumors may have global alterations in DNA methylation that could have pleiotropic effects on tumor–immune interactions. For instance, epi-MMRd tumors may have altered expression of immunomodulatory factors such as checkpoint molecules, which would thereby affect the state of the circulating immune response. It would be of great interest in future studies to longitudinally profile epi-MMRd tumors over the course of pembrolizumab therapy to investigate these and other potential mechanisms.

In other cancer types such as colorectal cancer, epi-MMRd tumors tend to develop at an older age compared with mut-MMRd tumors (75). Because the functionality of the immune system changes with aging (76), it is possible that some of the differential immune responses we observed in patients with epi-MMRd tumors are related to aging. However, we note that patient age was not significantly different between the epi-MMRd and mut-MMRd groups in our cohort, so the observed changes in immune system composition and gene expression are less likely to be due to age variations. That said, we cannot rule out that subsequent studies with larger patient cohorts may detect age differences between patients with epi-MMRd versus mut-MMRd endometrial tumors, which in turn may influence their response to anti–PD-1 immunotherapy (77).

An unexpected finding from our study was that pretreatment JAK1 mutations were not associated with primary resistance to anti–PD-1 ICB (Fig. 2B and C). At the same time, three patients had JAK1 mutations that were detected only after initiating pembrolizumab, a finding that is consistent with prior reports of acquired resistance to ICB through JAK1 alterations (24). We therefore hypothesize that the functional consequences of JAK1 mutations are partly determined by their temporal context, such that preexisting JAK1 mutations may reshape tumor–immune interactions in a manner distinct from when JAK1 mutations arise following ICB therapy. Interestingly, a recent study using an in vivo CRISPR screening approach demonstrated that pretreatment JAK–STAT mutations were associated with enhanced ICB sensitivity across multiple murine tumor models (78), in direct contrast to prior in vitro CRISPR screens that utilized coculture systems (25). Prior studies have also demonstrated in mouse models that loss of tumor-intrinsic IFNγ signaling can promote antitumor immunity through the enhanced activation of CD8⁺ T cells and NK cells (79–81). Mechanistically, tumors with IFNGR loss promote IFNγ production by CD8⁺ T cells, which in turn enhances the activation of cytotoxic NK cells (80, 81).

Although we had broadly classified the patients as nonresponders and responders based on standard RECIST criteria, we were also intrigued that some of the nonresponder patients did not experience disease progression for a year or more. These patients with long-term stabilized disease did not demonstrate substantial reductions in tumor size following pembrolizumab; nevertheless, we speculate that ICB treatment may have shifted the tumor–immune set point toward a new equilibrium, with the immune system restraining (but not effectively eliminating) the tumor. Future studies with larger patient cohorts should seek to interrogate the underlying tumor–immune interactions in patients with long-term stable disease.

Given the observational nature of the clinical trial design, the conclusions that can be drawn from the present study are largely correlative in nature. Nevertheless, we believe that a fruitful area for future research will be to investigate the underlying mechanisms driving the two distinct modes of antitumor immunity in patients with epi-MMRd versus mut-MMRd tumors. For instance, a deeper characterization of the tumor immune microenvironment in epi-MMRd and mut-MMRd tumors during ICB therapy would be highly informative. In addition to uncovering additional tumor-intrinsic features that differentiate the two classes of MMRd tumors, these data could reveal whether the distinguishing features of the circulating immune response we identified here are similarly observed in tumor-infiltrating immune cells.

Taken together, our findings fill an important gap in our understanding of the factors that influence antitumor immunity in the setting of a substantial neoantigen burden. By identifying the tumor-intrinsic and tumor-extrinsic factors that are associated with response to anti–PD-1 immunotherapy in MSI-H/MMRd endometrial cancers, our work illuminates new avenues for biomarker development and therapeutic intervention. Moving forward, it will be of interest to understand precisely how the molecular mechanism of MMRd can mold the antitumor immune response.

The trial was approved by the institutional review board of Yale University (HIC: 1605017712). Written informed consent was obtained from all participants. This study was conducted in accordance with the U.S. Common Rule and the Declaration of Helsinki. The study design and patient population have been previously reported in detail (11). In brief, NCT02899793 is an investigator-initiated phase II trial evaluating the safety and efficacy of pembrolizumab in patients with recurrent MMRd/MSI-H endometrial cancer, determined by IHC for MMR proteins and by microsatellite PCR. Patients were treated with pembrolizumab 200 mg intravenously every 3 weeks, until the occurrence of disease progression or adverse effects prohibiting continued therapy. The primary efficacy endpoint was objective response, as determined by RECIST v1.1 criteria. Secondary outcome measures were the duration of progression-free survival and overall survival. This study was partly supported by Merck, which supplied the study drug pembrolizumab but did not influence any aspect of the study design, data collection, analyses, or interpretation.

MSI-H/MMRd tumors with MLH1 promoter methylation were classified as epi-MMRd; tumors with somatic MMR gene mutations and without MLH1 promoter methylation were classified as mut-MMRd. In cases where a tumor exhibited both MLH1 promoter methylation and somatic MMR gene mutations, the methylation event was considered the primary driver, as epi-MMRd can drive the subsequent acquisition of MMR gene mutations due to heightened mutagenic potential (75, 82). As a confirmatory test in such cases, we leveraged the fact that tumors with epi-MMRd classically exhibit loss of MLH1 and PMS2 on IHC (83, 84). Thus, for tumors exhibiting both MLH1 promoter methylation and somatic MMR gene mutations, loss of MLH1 and PMS2 by IHC helped confirm the classification of epi-MMRd. Patient PEM06 was a special case in that the tumor exhibited complete loss of MSH6, MLH1, and PMS2 by IHC, along with MLH1 promoter methylation and biallelic loss-of-function MSH6 mutations. Given these features, we classified patient PEM06 as mut-MMRd despite the presence of MLH1 promoter methylation. Our reasoning was that if MLH1 promoter methylation was instead the initial driving event for MMRd in patient PEM06, and the MSH6 mutations were randomly acquired later due to the heightened mutagenic potential from epi-MMRd, we would expect MSH6 loss by IHC only in a subset of the tumor cells. The complete loss of MSH6 by IHC in PEM06 is therefore less consistent with epi-MMRd, so we inferred that MSH6 loss was the initial event driving MMRd (hence the mut-MMRd classification).

Genomic DNA was extracted from the tumor samples and captured on the NimbleGen 2.1M human exome array, then sequenced on an Illumina HiSeq 2000. WES data were aligned to the hg19 genome (human G1Kv37 reference with decoy sequences) using BWA-MEM (RRID:SCR_010910; ref. 85), with duplicate filtering performed by Picard (RRID:SCR_006525). The reads were then processed with GATK (RRID:SCR_001876; refs. 86, 87) according to the Best Practice workflow. Matched normal WES data were utilized as a reference to call somatic mutations using MuTect2 (RRID:SCR_001876), with DeTiN analysis to improve detection sensitivity (88). Predicted pathogenic/deleterious exonic mutations were defined as frameshift mutations, nonsense mutations, splicing site mutations, or missense mutations that were predicted to be pathogenic by at least two mutation assessors [out of MutationAssessor (RRID:SCR_005762; ref. 89), SIFT (RRID:SCR_012813; ref. 90), and CADD (RRID:SCR_018393; ref. 91)]. We refer to these mutations as “predicted pathogenic” throughout, because many of these mutations have not been experimentally validated. The predicted pathogenic/deleterious mutation calls were summarized and visualized using Maftools (92). To assess the association between specific gene mutations and ICB response, we used Fisher exact test on the 2 × 2 contingency matrix between mutation status and ICB response. A pseudocount was added when calculating the log odds ratio for ARID1A, as all patients with nonmutant ARID1A were responsive to therapy. HLA typing was performed on the matched normal WES profiles using HLA-HD (93). Using the high-confidence exonic mutation set above, we predicted potential MHC-I and MHC-II neoantigens in each sample using pVAC-Seq (94).

For a detailed analysis of JAK1 mutations, we transformed all exonic VAFs from MuTect2 into CCFs (39, 95). To do so, we first utilized Sequenza (RRID:SCR_016662; ref. 96) on the aligned .bam files to estimate tumor cellularity and to generate allele-specific copy-number segmentation profiles. We then calculated nonclustered CCFs using the standard formula (39), as well as clustered CCF estimates with PyClone-VI (97). For the purpose of visualization, when comparing the CCFs for the same JAK1 variant in pretreatment versus posttreatment samples, we assigned a CCF of 0 when the variant was not detected in a particular time point.

For the external validation of our findings regarding the relation between the mechanism of MMRd, TMB, and immunotherapy response, we reanalyzed two independent cohorts of patients with endometrial cancer. From the TCGA UCEC study (13), we categorized the patients into four groups based on MSI status and the mechanism of MMR deficiency. MSI/MMR status was based on the TCGA UCEC subtype annotations, as reported in a prior publication (14). Epigenetic MMRd was determined by analysis of methylation levels of the MLH1 promoter, regardless of the presence of MMR gene mutations. Specifically, we utilized the MEXPRESS methylation browser (98) to identify probes in the MLH1 promoter with a correlation coefficient ≤ −0.8 with MLH1 gene expression. The b-values for these 29 probes were then averaged within each sample to obtain a composite methylation b-value for the MLH1 promoter. We defined MLH1 promoter methylation with a cutoff of 0.5. Mutational MMRd was defined by the absence of MLH1 promoter methylation and the presence of mutations in any of the canonical MMR genes: MLH1, MSH2, MSH6, or PMS2. All other MSI-H tumors that were not classified as mut-MMRd or epi-MMRd were labeled as MMRd not otherwise specified (other). Finally, we utilized the nonsynonymous mutation rates and predicted neoantigen loads from a prior publication (14) to assess the relation between MMRd mechanism, TMB, and neoantigen load.

From the MSK-IMPACT cohort (15–17), we filtered for patients with endometrial cancer and categorized them into three groups based on MSI status and the presence of mutations in canonical MMR genes (MLH1, MSH2, MSH6, or PMS2). Of note, this data set did not have available annotations of MLH1 promoter methylation status. The three groups were thus (i) MMRp, (ii) mut-MMRd, and (iii) nonmutational MMRd (nonmut-MMRd). mut-MMRd patients were defined by the combination of MSI-H status and the presence of mutations in MLH1, MSH2, MSH6, or PMS2. In contrast, nonmut-MMRd patients were defined by the combination of MSI-H status and the absence of mutations in MLH1, MSH2, MSH6, or PMS2. We then utilized the authors’ annotations to assess the relation between MMRd mechanism, TMB, and ICB response.

The prevalence of JAK1 alterations across cancer types was queried using cBioPortal (99, 100), examining JAK1 mutations and homozygous deletions in the TCGA Pan-Cancer Atlas (36) and the AACR–GENIE cohort (38).

PBMCs were isolated by standard density gradient centrifugation with Ficoll-Paque (GE Healthcare) and used to construct scRNA-seq libraries. The Chromium Single-Cell 5′ and Chromium Single-Cell V(D)J Reagent Kit (10X Genomics) was used following the manufacturer's instructions.

scRNA/scTCR-seq reads were processed using Cell Ranger (RRID:SCR_017344; ref. 101) to generate both raw expression count matrices and TCR sequence information. Healthy donor PBMC scRNA-seq and scTCR-seq control data were obtained from prior publications: Pappalardo and colleagues (102) and Song and colleagues (103). Expression data were loaded into Seurat (RRID:SCR_016341; ref. 43) and prefiltered for cells with ≥400 counts, ≥150 unique expressed genes, and ≤20% mitochondrial gene expression. Cells were not filtered for complexity or doublet-related markers. The data were normalized with SCTransform (104), regressing out counts per cell, features per cell, and mitochondrial gene expression. The data sets were then integrated using reciprocal principal component analyisis (PCA) to identify integration anchors before performing cell clustering and UMAP dimensional reduction.

Clusters were annotated using manual analysis of marker gene expression in combination with automatic cell type identification using SingleR (105). Fourteen T-cell clusters were subsetted and reprocessed as a separate data set by rerunning PCA, UMAP, and clustering. After an initial round of analysis, two clusters from the T-cell data set were determined to be composed of cells with uniformly high mitochondrial reads. Further clusters were determined to be differentiated almost entirely by expression of a single TCR variable gene while encompassing many distinct T-cell subtypes. To address these issues, we regenerated the data set by prefiltering out the high mitochondrial read cells and excluding all TCR variables and joining genes from the PCA and clustering analyses.

To identify differentially enriched cell neighborhoods in the scRNA-seq data set, we utilized Milo (40). The single cells were first assigned to partially overlapping neighborhoods based on a k-nearest neighbor graph, with settings of k = 200 and d = 50 that were selected based on the distribution of neighborhood sizes, per the authors’ recommendations. After quantifying the number of cells in each neighborhood for each sample, we performed differential abundance testing with a negative binomial generalized linear model as implemented in edgeR (106). Multiple hypothesis correction was performed with a weighted FDR procedure that accounts for the spatial overlapping of neighborhoods. We chose a cutoff of spatial FDR (q) < 0.1 for identifying cell neighborhoods that were differentially enriched.

For characterizing the marker genes of differentially enriched cell neighborhoods, we calculated the average expression in each neighborhood for the 5,000 most variable genes in the data set. We then identified statistically significant marker genes by the two-tailed unpaired Mann–Whitney test with Benjamini–Hochberg multiple hypothesis correction, comparing enriched and nonenriched cell neighborhoods.

We identified high-confidence single-cell TCRαβ sequences using the Cell Ranger VDJ pipeline. The data were loaded into CellaRepertorium and restricted to TCRα and TCRβ chains with CDR3s that were high confidence (no N's), full length, ≥6 amino acids long, and productive. Data from bulk TCR-seq of TILs were loaded into Immunarch (10.5281/zenodo.6984421) and filtered for in-frame, coding CDR3s. Because the bulk data contained only TCRβ sequences, TIL-matching cells were identified in the PBMC scTCR-seq data by matching the scTCRβ CDR3s to the bulk TIL TCR-seq data from the same patient, similar to the approach utilized in a previous study (107). TCR analyses were therefore conducted on TCRβ sequences alone in order to maintain consistency and allow for the identification of tumor-matching TCRs. To determine the clonality of the TCR repertoire in each sample, we calculated the Simpson clonality index.

Conventional scRNA-seq differential expression analysis operates on the level of single cells, such as the Mann–Whitney test that is commonly used in the Seurat analysis package. However, this approach has a number of important limitations when analyzing large multisubject data sets: samples with larger cell numbers will dominate the statistical comparisons, and each individual cell will be counted as an entirely independent event, artificially inflating statistical power. To circumvent these issues, we used an alternative approach to differential expression analysis of scRNA-seq data. For each scRNA-seq library, we generated five random samples of each cell type and then aggregated the counts together. In effect, each cell type within a given library was condensed into five replicates. These “pseudobulk” replicates were then analyzed using DESeq2 (45), which has been demonstrated to have robust performance on scRNA-seq data despite being developed for bulk RNA-seq analysis (44, 108). DEGs were defined as those with adjusted P < 0.05. KEGG pathway analysis was performed with clusterProfiler (109) on the set of DEGs, using a hypergeometric test with Benjamini–Hochberg multiple hypothesis correction to evaluate statistical significance. Signed −log₁₀q-values were calculated by multiplying the −log₁₀q-values by the directionality of differential expression, such that pathways associated with downregulated DEGs would have a negative −log₁₀q-value. In the rare cases where a given pathway was associated with both upregulated and downregulated DEGs, we selected the directionality with the strongest statistical significance for visualization.

We obtained overall survival data (110) from the TCGA UCEC cohort, evaluating the relationship between survival and several factors previously associated with response to immunotherapy. These factors were previously computed across the TCGA data set in a prior study (14). We constructed a multivariable Cox proportional hazards model with the following variables, all coded as continuous variables: homologous recombination defects, nonsilent mutation rate, TGFβ response, IFNγ response, MSI score (MSIsensor, available on cBioPortal), CD8⁺ T-cell abundance, and activated NK cell abundance. The resulting hazard ratios were visualized as a forest plot.

To evaluate whether the gene signatures of CD16⁺ NK cells in epiR patients were associated with patient survival, we first compiled a final set of 111 genes that were concordantly upregulated or downregulated in epiR patients compared with NR and mutR patients. These 111 genes thus met several criteria: (i) they were called a DEG in epiR versus NR, both pretreatment and posttreatment, with the same directionality of differential expression; (ii) they were called a DEG in epiR versus mutR, both pretreatment and posttreatment, with the same directionality of differential expression; and (iii) the directionality of differential expression was shared between the epiR versus NR and epiR versus mutR comparisons. The 111 genes were utilized in a Lasso Cox regression to identify a core set of genes associated with survival. Of the five genes retained by the Lasso model, one gene was filtered out (PITPNC1), as high expression of this gene was associated with better survival, but this gene was downregulated in CD16⁺ NK cells. The four remaining genes were termed the epiR-NK4 gene set, comprised of CD63, PPIB, CEBPB, and LDOC1. In a multivariable Cox regression model, all four genes were independently associated with survival. To generate a summative score of the epiR-NK4 gene set, we performed PCA on the TCGA data set, based on the expression of the epiR-NK4 genes. We took the PC1 as the epiR-NK4 score and stratified patients into high versus low groups by the median score. Kaplan–Meier survival curves were then generated with the epiR-NK4 score classifications, and statistical significance was assessed by the log-rank test.

Raw and processed scRNA-seq/scTCR-seq data generated in this study are available in the Gene Expression Omnibus (GSE212217). Raw bulk TCR-seq data of TILs were generated by an external core facility (Adaptive Biotechnologies) and are unavailable, but the processed data are available in Zenodo (10.5281/zenodo.6985601) and the Gene Expression Omnibus (GSE212217). Raw exome sequencing data are unavailable due to the provisions of the trial sponsors; the processed exome variant calls and copy-number segmentation files are available in Zenodo (10.5281/zenodo.6985601) and GitHub (https://github.com/rdchow/endo_PD1). The predicted pathogenic/deleterious variant calls are also available in the supplementary tables. Custom scripts used for analysis and additional processed data files are available on GitHub (https://github.com/rdchow/endo_PD1). Previously published data sets analyzed in this study are described in the Methods.

R.D. Chow reports grants from the NIH/NCI and other support from the Paul and Daisy Soros Foundation during the conduct of the study. A. Iwasaki reports grants from Howard Hughes Medical Institute during the conduct of the study, as well as other support from RIGImmune and Xanadu outside the submitted work. E. Song reports grants from NIH Medical Scientist Training Program T32GM136651, NCI F30CA239444, and Paul and Daisy Soros during the conduct of the study. A.D. Santin reports grants from Puma, Immunomedics, Gilead, Synthon, Boehringer Ingelheim, and Genentech, grants and personal fees from Merck, Tesaro, and R-Pharm, and personal fees from Eisai outside the submitted work. No disclosures were reported by the other authors.

R.D. Chow: Conceptualization, data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. T. Michaels: Software, formal analysis, methodology, writing–review and editing. S. Bellone: Resources, data curation, investigation, writing–review and editing. T.M.P. Hartwich: Resources, data curation, writing–review and editing. E. Bonazzoli: Data curation, investigation, writing–review and editing. A. Iwasaki: Conceptualization, resources, supervision, funding acquisition, investigation, writing–review and editing. E. Song: Conceptualization, data curation, software, formal analysis, supervision, investigation, visualization, methodology, writing–original draft, writing–review and editing. A.D. Santin: Conceptualization, resources, supervision, funding acquisition, project administration, writing–review and editing.

R.D. Chow is supported by the NIH Medical Scientist Training Program (T32GM136651) and the NIH/NCI (F30CA250249). A. Iwasaki is an Investigator of the Howard Hughes Medical Institute. E. Song is supported by the NIH Medical Scientist Training Program (T32GM136651) and the NIH/NCI (F30CA239444). R.D. Chow and E. Song are also supported by the Paul and Daisy Soros Fellowship for New Americans. A.D. Santin is supported in part by Gilead Sciences, Inc., the NIH (U01CA176067-01A1), the NCI (P30CA16359), Stand Up To Cancer (Convergence Grant 2.0), the Discovery to Cure Foundation, and the Guido Berlucchi Foundation. The authors also thank Merck-US for its industry support.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.