Background

Uveal melanoma (UM) is a rare form of melanoma and the most common primary malignancy of the eye¹. It develops in the uvea of the eye, most often in the choroid and ciliary body, and more infrequently in the iris. Primary UM can be cured by brachytherapy or enucleation of the eye, but 50% of patients will develop metastasis². The most common route for metastatic spread is the liver for reasons that are poorly understood. Patients with metastatic disease have a median survival of less than a year ³; however, ongoing clinical studies suggest that some treatments could prolong survival ⁴. Monotherapy immune checkpoint inhibitors (ICI) are markedly less effective in patients with UM^5,6 than in those with metastatic cutaneous melanomas. However, combination treatments with PD-1/CTLA-4 inhibitors^7–9 or PD-1/HDAC inhibitors^10,11 have demonstrated longer overall survival than historic benchmark data ³. Tebentafusp, a T-cell engager suitable for patients with the HLA- A2 genotype^12,13, is the first therapy to show a prolonged overall survival in a phase 3 randomized trial, increasing median survival from 16.0 months to 21.7 months ¹⁴. Another development is locoregional treatments, where a recent phase 3 trial demonstrated that isolated hepatic perfusion (IHP) triples hepatic progression-free survival¹⁵ compared with the best alternative care and historic benchmark data. Although overall survival data are not yet mature, retrospective data suggest a benefit with liver-directed therapy¹⁶. However, invariably, all patients progressed on both tebentafusp and IHP; therefore, more research is needed.

Adoptive cell therapy (ACT) with TILs or CAR-T cells has not been extensively studied in UM. Pilot data from a clinical trial ¹⁷ showed that TILs can cause responses in patients, and we have previously showed that HER2 CAR-T cells can eradicate UM in xenografts ¹⁸. There is, however, a lack of robust ex vivo screening models and very few patient-derived xenograft (PDX) mouse models from metastases ¹⁹ to use in studies to improve ACT. Part of the lack of success for ACT may be attributable to a lack of tumor-reactive lymphocytes, possibly owing to lower mutation burden in UM compared to cutaneous melanoma ²⁰. However, both the clinical trial and previous analyses of TILs in UM liver metastases have suggested that at least some TILs are tumor reactive^17,21–23.

Defining tumor-reactive T lymphocytes (TRLs) among TILs would enable the identification of candidate biomarkers for selective expansion or TCR transgenics in cell therapy experiments in mouse models and cell therapy trials. To decipher how immunotherapy can be made more effective, the aim of this study was to identify and functionally validate TRLs in metastases of patients with UM. To this end, we used paired single-cell RNA (scRNA) and TCR sequencing, as well as functional experiments using patient-derived xenograft (PDX) models and autologous TIL cultures.

Results

We previously collected metastases of UM from patients in routine clinical care and two clinical trials. These were the SCANDIUM phase 3 randomized trial comparing IHP to best alternative care²⁴ and the PEMDAC phase 2 single-arm trial where patients received a combination of the PD-1 inhibitor pembrolizumab and the HDAC inhibitor entinostat ¹¹. Owing to the limited amount available of these biopsies, we prioritized one part of the biopsy for DNA/RNA preparation to sequence²², one part for transplantation into immunocompromised mice to generate PDX mouse models, one part for generation of TILs, one part for cryopreservation of finely minced tumor (flow cytometry and single-cell sequencing) and one part for formalin-fixed paraffin embedded blocks (FFPE) (Fig. 1a).

Figure 1.

Single-cell sequencing defines an atlas of T cell states in UM metastases

We obtained biopsies suitable for sequencing from 14 patients, two of which were subcutaneous and 12 were liver metastases. To identify and profile the T cells present in these samples, we performed paired scRNA and TCR sequencing, integrated the samples, and annotated the cell types. This has resulted in the largest single-cell atlas of metastatic UM to date, complementing prior efforts, including our own study of TILs from UM metastases ^21,22,25,26.

Uniform Manifold Approximation and Projection (UMAP) dimensionality-reduced expression profiles revealed clusters containing melanoma, CD8⁺ or CD4⁺ T cells, NK cells, monocytes, and small clusters of B cells, endothelial cells, and other minor cell populations (Fig. 1b-c, Supplementary Fig. 1a-f). Clusters were generally well-mixed with respect to contributions from different samples, suggesting few patient-specific or batch effects after data integration (Supplementary Fig. 1g). The identity of cells labeled as UM was further confirmed by inferring copy number changes. This revealed characteristic aberrations seen in metastatic UM, including monosomy of chromosome 3 and gain of 8q, as well as subclonal variation within patients (Supplementary Fig. 1h).

Further sub-clustering of CD8⁺ T cells showed the presence of 12 different groups (Fig. 2a-b and Supplementary Fig. 1i-k), all characterized by the differential expression of immunological marker genes (Fig. 2c-d and Supplementary Table 1). The clusters can be divided into two axes. One of these is characterized by expression of genes associated with late activation and exhaustion, such as PDCD1, HAVCR2, TNFRSF9 and HLA- DRA. Additional expression of ITGAE (CD103) suggests that these clusters contain tissue-resident memory cells²⁷. Multiplex immunofluorescence showed that PD-1- and ICOS-expressing T cells were in close proximity to the tumor cells (Fig. 2e). The other axis expressed genes associated with naïve/memory-like or early activated phenotypes, such as IL7R, TCF7, CCR7 and NR4A1 (Fig. 2d).

Figure 2.

Individual clusters within these axes mostly correspond to intermediate phenotypes that progress from activation to exhaustion. However, a few seem to represent distinct states. Two of these were within the IL7⁺ axis. Of these, one expressed KLRB1 and NCR3 and the other expressed NK-associated markers, such as FCGR3A, FCRL6 and NKG7 (Fig. 2c-d). Although both clusters expressed genes related to cytotoxicity, only the latter expressed GZMB and GZMH. Other markers associated with this group are CX3CR1 and FGFBP2. Previous studies labeled these cells as cytotoxic NK-like T cells^28,29.

In summary, this map reveals the cell types present in liver and subcutaneous metastases of UM and divides infiltrating T cells into two major phenotypic groups, characterized by the mutually exclusive expression of HAVCR2 and IL7R, but also reveals distinct cytotoxic subsets expressing NK-associated markers.

TILs from UM metastases can kill their cognate tumors in PDX models

Given the presence of activated T cells in tumors, we investigated whether these could eradicate metastatic UM in ACT experiments using PDX models. These models are useful personalized models for studying tumor biology and more recently, tumor immunology³⁰. However, the biopsies we received were generally very small, and the initial take rate using the same single-cell dispersion method used for cutaneous melanoma metastases³⁰ was unsatisfactory. Therefore, we transplanted small tumor pieces or used splenic injections, which led to successful PDX generation for some patients. Ten of these are presented in this study. All the PDX tumors exhibited confirmatory markers of melanoma, such as SOX10, MELAN-A (MART1), HMB45 or S100B. Sanger sequencing further revealed mutually exclusive mutations in the UM oncogenes GNAQ, GNA11, CYSLTR2 and PLCB4, as found in biopsies²². (Fig. 3a and Supplementary Fig. 2). The UM tumors grew at a slower rate compared to our experience with cutaneous melanoma, which has a take rate of 95% and where a model is usually established within 3-12 weeks (Fig. 3b)³¹.

Figure 3.

We have previously developed a method to study ACT in PDX models, termed PDXv2^18,32. Survival, expansion, and activity of TILs and CAR-T cells were improved in this model using human IL2-expressing transgenic NOD SCID IL2 receptor gamma knockout mice (hIL2-NOG)^18,32–35. Tumor eradication by TILs correlates with immunotherapy response in patients with melanoma^32,36. We conducted a PDXv2 experiment using PDX models from two patients with UM. The injection of autologous TILs into hIL2-NOG mice carrying UM1 resulted in tumor rejection, whereas UM22 tumors could not be eradicated (Fig. 3c-d). This constitutes proof of concept that TILs from UM metastases can kill autologous tumors in a PDXv2 model.

TILs contain subsets capable of ex vivo expansion and killing in 3D tumor sphere cultures

To gain further insight into the capacity of T cells in UM biopsies to kill tumors, we next investigated the tumor reactivity of ex vivo expanded TILs on autologous UM cells. Since most UM PDX models did not grow fast enough to routinely conduct ACT experiments, we used the tumors grown in the mice to generate 3D tumor sphere cultures ^37,38 of high viability for TIL experiments^39,40. We generated tumor sphere cultures from ten PDX models and expanded their autologous TILs (Fig. 3e). When tumor spheres and TILs were co-cultured, TILs from patients UM1 and UM9 were activated the most (surface 4-1BB⁺) and also degranulated (surface CD107a⁺; Fig. 3f-g). These TILs also bound to the tumor spheres, as assessed by labeling tumor spheres and TILs with different vital dyes and then imaging with a real-time cell culture imager (Supplementary Fig. 3a-b). The UM22 TILs were hardly able to get activated by their autologous tumor spheres (Fig. 3b), which likely explains their lack of activity in the PDX ACT experiment (Fig. 3d). However, since UM22 is capable of presenting MART1 on HLA-A2, and since we previously identified a minor population of MART1-reactive TILs in this biopsy⁴¹, we decided to expand and enrich this subset using MART1 peptide HLA-A2 dextramers. When using these cells, we were able to achieve potent activation and degranulation (Fig. 3f-g). This confirms that T cells extracted from biopsies can contain subsets capable of recognizing and killing tumors, which can be expanded ex vivo.

Determining the phenotypes of TRLs in biopsies using scRNA and TCR sequencing

To gain clues on how to identify reactive T cells from biopsies through their transcriptomic phenotypes, we co- cultured UM1 and UM9 TILs with their autologous tumors and sorted 4-1BB activated T cells by cell sorting. We prepared RNA from these cells, and as a control, from MART1-specific allogenic UM46 TILs and sequenced their TCRs. Then, using the single-cell atlas built from biopsies, we mapped those TCRs back to T cell subsets from their original tumors. This showed that those TRLs predominantly emanated from exhausted and late-activated cells (Fig. 4a-b). UM9 also contained a sizeable amount of TRLs in the NK-like/cytolytic cluster (Fig. 4b). The lack of potential TRL clones identified in the TCF7⁺ axis of the CD8⁺ T cell atlas from biopsies is also compatible with a previous report showing that TCF7⁺ T cells in melanoma tend to be nonreactive bystanders²⁸.

Figure 4.

To evaluate which T cells in biopsies were reflective of those capable of generating TIL cultures, we further mapped our previously published TCR sequencing data of UM TILs to the biopsy data (Fig. 4c) to compare them to the clones matching experimentally identified TRLs in the biopsy (Fig. 4b-d, Supplementary Fig. 4a-b and Supplementary Table 2). Cells from the TIL dataset tended to correspond to biopsy T cells with exhausted and late-activated profiles, similar to the experimentally identified TRLs. This might be expected since both have undergone expansion, which might favor T cells in certain phenotypic states. However, TRLs were significantly over- and underrepresented in some phenotypic clusters compared to TILs in general. Overrepresented clusters included proliferative exhausted and IL7R⁺ exhausted T-cells. The underrepresented ones were activated (non- exhausted), naïve/memory-like, and FCGR3A⁺ NK-like T cells. The distribution of TRLs across the clusters was also different for cells sorted by 4-1BB-positivity versus MART1-recognition. MART1-reactive cells were more likely to be found in exhausted and late-activated clusters than 4-1BB⁺ cells, which tended to map to earlier activation states as well to the NK-like cluster (Fig. 4e). The latter may be partly due to their prevalence in UM9 (Supplementary Fig. 4c).

Mapping TCRs from enriched 4-1BB⁺ cells to TILs in corresponding patients from previously published TCR sequencing data revealed that a large fraction of TILs in these cultures was among the identified TRL clones (Fig. 4f). In 6 of the 8 previously published TIL cultures, we were able to recover and sequence the corresponding cryopreserved biopsy of that patient. Five of the six cultures had many TCRs that overlapped between TILs and biopsies (Fig. 4g and Supplementary Fig. 4a-b). UM3 did not have any matching TCRs between the biopsies and expanded TIL cultures. No overlap was observed between samples from different patients with the 4- 1BB⁺/MART1⁺-selected TCRs (Supplementary Fig. 4d). This is reasonable given the different HLA genotypes of these patients. In general, very few clonotypes were shared between samples from different patients (Supplementary Fig. 4e). In summary, experimentally identified TRLs mapped back to biopsies were overrepresented in proliferative exhausted and IL7R⁺ exhausted T cells, and these T cells could also be identified in sequenced cultures of expanded TILs.

Identifying potential antigens recognized by TRLs

To understand potential antigens recognized by a given TCR, we searched for matches in public databases of known TCRβ chains, based on CDR3β motifs, using TCRMatch⁴². The results confirmed that a number of TCRs matched known MART1-reactive clonotypes but also revealed matching TCRs that might recognize other melanoma antigens, such as PMEL and MAGEA1 (Supplementary Fig. 4f-g, h). Overall, a large number of TCRs also matched clonotypes recognizing viral antigens such as influenza, EBV, and human coronaviruses. However, an even greater number lacked known matches, suggesting they might recognize tumor-specific antigens (neoantigens).

However, public databases are currently limited. As an alternative approach to nominate tumor- reactive clonotypes, we used an algorithm designed to predict clusters of TCRs that are sufficiently similar and HLA-restricted to potentially recognize the same antigens. GLIPH2⁴³ nominated a few high-confidence clusters of similar TCR sequences based on the biopsy scRNA-seq data (Fig. 4h and Supplementary Table 3). To further test the hypothesis that each cluster recognizes the same antigen, we mapped the CDR3β motifs of these sequences to the same public databases. Interestingly, the highest-ranked cluster was overrepresented for sequences matching MART1, and this cluster also exhibited significant restriction towards HLA-A*02:01, on which MART1 is known to be presented. The remaining high-confidence clusters mainly matched antigens of viral origin or lacked matches to known antigens (Fig. 4h).

Similar to the experimentally determined MART1-reactive T cells, the bioinformatically determined cluster of potential MART1-recognizing clones mostly belonged to the exhausted CD8⁺ T cell cluster (Fig. 4i, Supplementary Fig. 4i). This supports the finding that TRLs found in biopsies may have predominantly exhausted profiles. In addition, a high-confidence group of clonotypes co-localized to the NK-like cluster was nominated by GLIPH2, although no sequences from this cluster matched the public TCR-antigen pairs (Fig. 4h, Supplementary Fig. 4i). Given that the cells from this cluster were experimentally identified as 4-1BB⁺ upon tumor co-culture (Fig. 4b), one might speculate that they could recognize unknown tumor antigens. Taken together, these findings suggest that both bioinformatically inferred TRLs from biopsies and ex vivo expanded TRLs tend to have exhausted phenotypes. Among these, T cells recognize known melanoma antigens such as MART1, PMEL, and MAGEA1, as well as unknown sequences that might be neoantigens.

TRLs home to the tumor of PDX mice bearing liver metastases of UM

TIL cultures contain both TRLs, virus-specific T cells, and bystander T cells that grow during IL2-containing culture conditions. We had previously found that TILs can home to subcutaneously growing skin melanoma metastases ^18,32 but we were interested in determining whether they could identify and home to UM growing in the liver. The PDX models used to generate sphere cultures were prepared by subcutaneous transplantation. To instead generate a liver metastasis model, we performed splenic injection of UM22 into NOG mice⁴⁴. Tumor cells forming liver metastases were further transplanted via tail vein injection into hIL2-NOG mice. After verifying tumors by ultrasound, autologous bulk or MART1-selected UM22 TILs were intravenously injected for TIL therapy studies (Fig. 5a). The mice were harvested three weeks later, and tumors recovered from the liver were analyzed by flow cytometry, IHC, single-cell RNA and TCR sequencing. As expected from Fig. 3f-g, MART1-selected TILs had a higher degree of activation than unselected bulk UM22 TILs (Fig. 5b), despite similar numbers of TILs present in the livers. There were no visible differences in the size of the liver tumors between mice receiving unselected or MART1-selected UM22 TILs. This may be due to T-cell exclusion, since most TILs colonized the rim of the tumor (Fig. 5c). More MART1-specific TILs than unselected TILs were able to get activated, as assessed by the activation markers 4-1BB and CD69 (Fig. 5b and d-e). This shows that T cells recognizing a tumor antigen can successfully migrate to liver metastases of UM in mice.

Figure 5.

Determining the phenotypes of TRLs in PDX liver metastases using scRNA and TCR sequencing

Next, we profiled the livers and spleens of mice injected with unselected bulk or MART1-selected UM22 TILs by scRNA and TCR sequencing (Fig. 5f, Supplementary Fig. 5a). After bioinformatically separating human- from mouse-derived cells and integrating the data, clearly distinguishable clusters emerged containing UM cells, all characterized by the expression of the differentiation genes TYR, MLANA and PMEL. Their designation as UM cells was further confirmed by inference of their copy number profiles (Supplementary Fig. 5b). As expected, all UM cells resided in the liver (Supplementary Fig. 5a), despite the mice being injected with UM cells via the tail vein. The TILs discovered in these samples constituted three clusters: C1, C2, and a proliferative cluster with a mixed set of cells from the two main clusters (Fig. 5g). Compared with CD8⁺ T cell clusters from the biopsies, the clusters found in mice did not separate with respect to the same set of markers. Rather, all clusters expressed markers compatible with late-activated or exhausted cytotoxic cells to some extent (Supplementary Fig. 5c). These TILs were detectable in both the liver and the spleen, although there was a notable correspondence between, on one hand, liver TILs and C1, and on the other, spleen TILs and C2 (Fig. 5g-i).

Investigating the TCR clonotypes represented in the four experimental conditions revealed several clones unique to or overrepresented in the liver biopsies of mice with MART1-selected T cells, as compared to unselected TILs or spleen samples (Fig. 5j). Given that these clones have been sorted to enrich for reactivity towards an antigen known to be presented by UM22 cells, and that they have successfully migrated to the liver tumor, this subset is likely to be strongly enriched for TRLs. Among these, some were also found in previous scRNA-seq data from expanded UM22 TILs (Supplementary Fig. 5d-i).

Differential expression analysis between this TRL subset and unselected TILs found in the livers of PDX mice revealed statistically significant genes that may be potential markers of the tumor-reactive phenotype (Fig. 5k, Supplementary Fig. 5j, Supplementary Table 4). This included the activation marker TNFRFS9 (4-1BB) as well as KIR2DL4, which has been found to be expressed in certain exhausted subsets of CD8⁺ T cells (Fig. 2c)²⁸. Interestingly, this subset also bears a striking resemblance to a recently described group of TIM-3⁺IL7R^- tumor- reactive tissue-resident memory cells in lung cancer, further characterized by the specific expression of KRT86, DUSP4 and LAYN (Supplementary Fig. 6)⁴⁵. This set was enriched in patients responding to anti-PD1 therapy⁴⁵. Cells expressing this marker combination were also present in biopsy samples, predominantly in exhausted clusters (Fig. 5l-m), as well as in TIL cultures from UM22 (Supplementary Fig. 6b). This might suggest avenues for enriching TRLs from biopsies for use with ex vivo expansion and ACT based on their gene expression profiles.

CD8 cells diminish after tumor eradication in PDX models

Although the UM22 PDX model demonstrated signs of TIL homing and activation, TILs were unable to cure the mice from the tumor, neither when grown subcutaneously (Fig. 3d) nor in the liver (Fig. 5c). We therefore developed a liver metastasis model of UM1, which responds to TILs in a subcutaneous PDX model (Fig. 3d). Injection of TILs into hIL2-NOG mice bearing UM1 liver metastases resulted in a complete clearance or tumors as assessed macroscopically and by ultrasound imaging. To characterize the cellular composition of T cells in the livers and spleens of the PDX mice treated with TILs, we performed single-cell RNA and TCR sequencing. Clustering analysis identified 15 clusters of human cells that were annotated by known marker genes defined by differential expression analysis and supported by literature (Fig. 6b-c, Supplementary Fig. 6c-e, Supplementary Table 5 and 6). Notably, and at variance with the UM22 model, none of the clusters contained any melanoma cells, confirming that the TILs had eradicated the tumor.

Figure 6.

We compared number of cells coming from liver and spleen samples. A total of 20,418 cells were found in liver and 18,208 from spleen (Supplementary Fig. 7a-b). Strikingly, even though the TILs injected were primarily CD8 positive, the T cells remaining in liver and spleen after tumor eradication were primarily CD4+ positive cells (Supplementary Fig. 7b-c). Most CD8 cells resided in Cluster 11 with fewer cells also in clusters 8 and 12 due to expression of cell cycle genes. Re-clustering of cluster 11 resulted in two subpopulations (11.0 and 11.1) annotated as CD8+ cytotoxic expressing killer cell lectin like receptor G1 (KLRG1) and CD8+ cytotoxic expressing higher levels of CD8A and CD8B. (Supplementary Fig. 7d). These cells also expressed higher levels of CTLA4, TIGIT and TNFRSF9 (encoding activation marker CD137/4-1BB) and the ψ8Tcell markers TRGC2/TRGV9. Cluster 8 expressed high levels of Ki-67 (MKI67) assembly factor for spindle microtubules (ASPM), stathmin 1 (STMN1), DNA topoisomerase II alpha (TOP2A), TPX2 microtubule nucleation factor (TPX2) corresponding to the proliferative cellular state. Cluster 12 expressed minichromosome maintenance complex component 2, 6, 7 (MCM2, MCM6, MCM5, MCM7), GINS complex subunit 2 (GINS2), helicase, lymphoid specific (HELLS) that are associated with the proliferation and cell cycle state.

Next, we investigated the CD4+ cells and found that clusters 1, 3, 5, 9 expressed programmed cell death 1 (PD-1, also, known as PDCD1), lymphocyte activation gene 3 (LAG3), cytotoxic T lymphocyte-associated protein 4 (CTLA-4), T cell immunoreceptor with immunoglobulin and ITIM domains (TIGIT), thymocyte selection associated high mobility group box (TOX) and eomesodermin (EOMES) marker genes that are associated with a dysfunctional or exhausted state and was annotated as a CD4+ exhaustion TOX+, CD4+ exhaustion PTCD1+, CD4+ exhaustion LAG3+ and CD4+ exhaustion TIGIT+ respectively. Three adjacent clusters 0, 2, 10 expressed interleukin 7 receptor (IL7R), rhotekin 2 (RTKN2), C-C motif chemokine receptor 3 (CCR3) and protein tyrosine phosphatase receptor type C (PTPRC, also known as CD45) marker genes corresponding to a memory T cells phenotype and were annotated as a CD4+ memory RTKN2+, CD4+ memory CCR3+ and CD4+ memory PTPRC+ respectively. More liver cells contributed to the CD4+ exhaustion cluster whereas the CD4+ memory cluster was primarily made up of spleen cells (Supplementary Fig. 7b). Finally, three adjacent clusters 7, 6, 4 were associated with a CD4+ Th17 phenotype that play important role in the inflammatory response expressing the transcription factor retinoic acid orphan receptor (ROR)γt (RORC), chemokine CC receptor 6 (CCR6), a killer cell lectin like receptor B1 (KLRB1, also known as CD161), interleukin 4 induced 1 (IL4I1), cyclin dependent kinase inhibitor 1A (CDKN1A), aquaporin 3 (AQP3) genes (Supplementary Fig. 7e) and annotated as a CD4+ Th17 CD161+ (KLRB1), CD4+ Th17 PTPN13+ and CD4+ Th17 PTPRM+ respectively.

Finally, we investigated the TCR clonotype abundance and visualized top ten the most abundant clonotypes (Fig. 7a). Clonotypes 302, 11155 and 21797 were the most abundant, and as expected, they were present in CD4 + cells. Clonotypes 302 mapped to the CD4+ exhaustion site, 111555 to the CD4+ memory and 21797 to the CD4+ Th17 cluster (Fig. 7b). To investigate if there were any differences in gene expression between these clonally expanded T cells in liver versus spleen, we performed pseudobulk analysis. Genes such as FASLG, ICOS, TIGIT, IL2RA that were higher in T cells residing in liver are known to be expressed during activation and exhaustion. (Fig. 7c). We also mapped where TRLs are residing by extracting the top ten T cell Receptor Beta (TRB) sequences obtained from the CD137 sorted UM1 TILs (Fig.3f-g and Fig.4). We identified a few of these in the PDX data and mapped them on the UMAP (Fig. 7d). TRLs were residing on the CD4+ Th17 and proliferation/cell cycle clusters.

Figure 7.

Discussion

UM continues to be a condition with unmet medical need, despite recent advances in locoregional therapy and immunotherapy⁴⁶. However, the fact that a T cell engager is prolonging survival in patients with UM having an HLA-A2 genotype¹⁴ and that immune checkpoint inhibitors in various combination therapies on occasion can yield clinical benefit ^7,8,11 suggests that if we learn more about the tumor immunology of UM, then we can devise new immunotherapies. Here, we generated datasets and models that can be used not only to characterize tumor immunity, but also to functionally test hypotheses and targets in autologous culture and in vivo systems.

ACT with TILs has also been tested in patients with metastatic UM¹⁷. Although there were responses in a few patients, the data clearly show room for improvement. The HAITILS trial (NCT04812470) addresses whether hepatic arterial infusion results in local accumulation of TILs in liver metastases. However, this will not circumvent the low mutational burden of UM (other than iris melanoma²²) and the potentially low number of TRLs in the cell product (that we see in this study). A greater understanding of TRLs in patients’ metastases is therefore of interest to learn how to enhance the activity of TIL products and how to expand, for instance, neoantigen- or melanoma-associated antigen-specific TRLs^47,48. Here, we aimed to identify tumor-reactive lymphocytes in terms of marker characteristics and functional assessments of TILs. These analyses suggest that some of the TILs in our cultures that became activated (4-1BB positive) when binding tumor cells were already TRLs (4-1BB or PD-1 positive) when they resided in the tumor. However, they were never able to clear the tumor, and most likely became exhausted or dysfunctional in the microenvironment.

Presence of TRL-specific markers adds validity to a strategy whereby enrichment of TRLs could be done directly from the tumor by sorting for those markers. Indeed, previous studies have employed activation, exhaustion, or tissue-resident markers such as 4-1BB⁴⁹, PD-1^50–52, CD39,^52,53 and CD103,^54,55 or the specificity of a TCR for a melanoma-associated or melanoma-specific antigen^56,57. However, a risk with this strategy is that some TRLs are missed. Most notably, we revealed TRLs in clusters of T cells negative for activation and exhaustion markers. These cells are interesting since they express high levels of perforin and granzyme B, suggesting that they are cytotoxic. It is tempting to speculate that these cells are effector or effector memory cytotoxic T cells since they express RNA encoding CD16, FCRL6, and CX3CR1^58–60. In the UM9 sample, most of the T cells in the biopsy that mapped to TRLs identified in the tumor sphere assay resided in this cluster.

PDX models have been developed for UM, but to our knowledge, this is the first time they have been used to study immunotherapy using autologous TILs and single-cell sequencing. Although one model responded to TIL therapy, the other model did not. In the UM22 liver metastasis model, we found evidence that T cells were pushed out of the tumor. This may partly explain the failure of immune checkpoint inhibitors in many patients ⁶¹ and why ACT with TILs does not work for the majority of patients with metastatic UM¹⁷. Nevertheless, we profiled TRLs that preferentially migrated to liver metastases in mice, as opposed to non-reactive T cells and those that migrate to the spleen. We identified marker genes and phenotypic states that might be used to identify relevant subsets present in patient biopsies and TIL cultures. These T cells appeared to resemble a recently described IL7R^-HAVCR2⁺KRT86⁺DUSP4⁺LAYN⁺ subset, which has been associated with response to anti- PD1 therapy in lung cancer⁴⁵. They also tended to express ITGAE (CD103), suggesting that they may be tissue- resident memory cells. This subset was identified in both PDX liver metastases, cultured TILs, and patient biopsies. The phenotypes of cells with this marker combination were predominantly exhaustion-like in the PDX samples as well as in biopsies. This is compatible the results from the profiling of biopsy CD8⁺ T cells, where mapping of experimentally identified MART1-reactive T cells back to biopsies and computational prediction of tumor antigen-recognizing T cells suggested an over-representation of TRL candidates in exhausted and late- activated gene expression-based clusters. Whether cells expressing these canonical exhaustion-related markers are truly dysfunctional or capable of eliciting an antitumor response remains to be explored⁴⁵.

In the UM1 liver metastasis model, the tumor was eradicated, demonstrating that TILs can cause a complete response also of liver metastases. The single-cell data obtained confirmed the loss of melanoma cells but it also resulted in a loss of CD8+ T cells. It is well known in immunology that the fate of CD8 T cells after an acute immunological reaction is either death of the effector cells or development of memory. A strength of the PDX experiment is that we could model. But since we harvested the livers after tumor eradication, a limitation is that we could not investigate the expression of genes in T cells that are actively combating the tumor. However, since the CD4+ T cells were present in high numbers, we could compute differences between the T cells of the same clonotype residing either in liver or in spleen. This analysis did indeed suggest that some level of activation/exhaustion remained even after tumor clearance. Moreover, some of the TCRs from the coculture experiment of tumorspheres and TILs designed to reveal TRLs, were indeed identified in vivo among the CD8 and the proliferative cells. We are therefore confident that PDX models such as UM22 and UM1, will be useful in future studies where T cells are sorted for potential markers prior to expansion or infusion, where the T cells are genetically engineered or grown in cocktails of factors to enhance potency and stemness. We also propose the model to be appropriate for further humanization to unravel the interplay between tumor cells, liver cells and stroma cells, myeloid cells and other immune cells.

When mapping TCRs to TCR databases, we noted not only similarities to known melanoma-associated antigens, but also to viruses. Importantly, TCRs that recognize different antigens can share CDR3β chains. The TCR alpha chain is even more promiscuous but not always reported in these databases, which makes it difficult to find accurate TCR pair matches. Therefore, until current databases have been expanded and improved, we cannot be certain about the exact antigen matches of some of these receptors. Nevertheless, an increasing body of literature suggests that bystander cells in tumors are often antigen-experienced, can be directed against bacteria or viruses, and even participate in tumor control, for instance during immune checkpoint inhibitor therapy^62–64. Whether the presence of foreign antigen-specific T cells is beneficial for cellular immunotherapy remains to be evaluated.

In summary, we profiled the tumor microenvironment of metastatic uveal melanoma and identified the phenotypes of tumor-reactive T cell subsets. The results of this study provide new avenues for targeted expansion of T cells to improve cellular and immune checkpoint immunotherapy. Importantly however, to potentiate TIL therapies in UM, immune cold tumors likely need to be turned hot, i.e., to create an inflammatory environment promoting T cell infiltration. Simultaneously, lymphodepletion is needed to avoid existing T cells sequestering IL2 needed for the TIL infusion product. Tumor-targeting treatments such as isolated or percutaneous hepatic perfusion (IHP and PHP) may fulfill these needs and could act as adjuvant therapy for both immune checkpoint inhibitor treatment and TIL therapy. As melphalan triggers immunogenic death, influx of new T cells into cold tumors may occur, partly explaining the very high response rate of these therapies¹⁵. Ongoing clinical trials are combining locoregional therapy with immune checkpoint inhibitors (NCT04463368) to capitalize on this concept. However, data from these trials have not yet been published. Combination of liver- directed therapy and TIL therapy is also warranted in future animal experiments and clinical trials. We propose that PDX models are optimal to capture the intra-patient heterogeneity of cancer and that multiple mechanisms in tumor immunology and treatment can be modelled at single-cell resolution.

Methods

Patient samples

According to the ethical approval of the Regional Ethical Review Board (#289-12, #44-18 and #144-13), the patients provided oral and written information with signed informed consent agreements. Biopsies were extracted from either subcutaneous or liver metastases. Tumor biopsies were divided into pieces that were snap-frozen, minced, embedded in formalin blocks, or used for tumor-infiltrating lymphocyte cultures and cryopreservation. Snap-frozen tumor fragments were homogenized using a Bullet Blender (Next Advance, NY, USA). DNA and RNA were extracted using the AllPrep DNA/RNA kit (Qiagen, Germany).

Immunostainings of sections

Formalin-fixed, paraffin-embedded (FFPE) patient samples were evaluated using H&E staining and IHC. IHC was performed using an autostainer (Autostainer Link 48, Dako) with primary anti-human antibodies against Melan- A (clone A103; Dako), S100 (IR504; Dako), PRAME (E7I1B, CST), SOX10 (E6B6I, CST), HMB45 (M0634; Dako), CD3 (D4V8L, CST), and CD137 (E6Z7F, CST). HRP DAB or HRP Magenta (Dako) was used to stain the proteins of interest, and counterstaining was performed using hematoxylin.

We also performed multiplex immunofluorescence using The PhenoCycler-Fusion system (Akoya Biosciences) at the Spatial Proteomics Facility at SciLifeLab, Stockholm. The data for CD3, ICOS, and PD-1 primary antibodies are shown here; however, the full dataset will be published elsewhere.

PDX models

Animal experiments were conducted in conformity with E.U. directive 2010/63 (Regional Animal Ethics Committee of Gothenburg approvals no. 2014-36 and 2018-1183). PDX models were generated by transplanting small pieces into the flank (subcutaneous) of immunocompromised, non-obese, severe combined immune deficient interleukin-2 chain receptor γ knockout mice (NOG mice) or hIL2-NOG mice (Taconic, Denmark). Cryopreserved single cells from biopsies were thawed and transplanted by splenic injection to create liver metastases in NOG mice or hIL2-NOG mice (Taconic, Denmark). All PDX tumors were analyzed using IHC, and their identity was verified by staining patient biopsies using clinically graded antibodies. For subcutaneous implants, tumor growth was monitored twice a week following TILS treatment. After establishment, the liver metastasis models formed liver metastases equally well by splenic injection and tail vein injection. The mice were evaluated using ultrasound scans before TIL treatment and grossly assessed at the time of sacrifice.

Generation of TILs and MART1-specific T cells

The metastatic fragments were placed in culture medium (90% RPMI 1640 (Invitrogen), 10% heat inactivated Human AB serum (HS, Sigma-Aldrich), 6000 IU/ml recombinant human IL2 (Peprotech), penicillin and streptomycin (Invitrogen) to generate young TILS (yTILs) as previously described ³². MART1 specific yTILs were identified as previously described ²² and sorted using FACSAria III (BD Biosciences). yTILs were expanded using a standard REP (Rapid Expansion Protocol), constituting irradiated (40 Gy) allogeneic feeder cells, CD3 antibody (clone: OKT3) (30 ng/ml) (Miltenyi), medium (AIM-V, Invitrogen) supplemented with 10% HS and 6000 IU/ml IL2 and monitored for 14 days with media changes ³². After completion of the REP, MART1 specificity was confirmed using MART1-specific dextramers (Immudex, Copenhagen, Denmark). REP and REP MART1 TILs were used for downstream, cell line, spheroid, and in vivo experiments.

Spheroid assay

PDX models were harvested, and tumor tissue was digested using a human tumor dissociation kit (Miltenyi) according to the manufacturer’s instructions in conjunction with the gentleMACS Octo tissue dissociator (Miltenyi). This single-cell suspension was used to prepare spheroids using an ultra-low attachment (ULA), 96- well, clear, round-bottomed plates (Corning,NY), and GBM-MG serum-free media (CLS cell lines service, Germany) for 72-96 hours in a 37°C incubator. To monitor treatment responses, cells were labeled with CellVue™ Jade Cell Labeling Kit (Thermo Fisher Scientific) and/or CellVue Claret Far Red Fluorescent Cell Linker Mini Kit (Sigma) and used for co-culture experiments. The 3D images were captured using IncuCyte® (Essen BioScience) +/- TILs treatment, followed by imaging and analysis at 72 h with IncuCyte® Live Cell Analysis. Tumor spheres stained with a far-red dye were measured for the red object area/µm² mask with >100 μm. Flow cytometry was used to confirm the positively stained immune cells (Jade-FITC channel) and cancer cells (far red-APC channel).

Co-culture and TCR sequencing

A human UM cell line derived from patient UM22 ²² was grown in complete medium (RPMI-1640 supplemented with 10% FBS, glutamine, and gentamycin) and cultured at 37°C with 5% CO2. The cell line or spheroids were analyzed with or without TILs, and the co-culture was analyzed using human CD3 (FITC and HIT3a, BioLegend), CD8 (APC and SK1, BioLegend), CD107a (APC and H4A3, BD Biosciences), CD137 (PE and 4B4-1, BioLegend), active caspase-3 (clone C92-605; BD Biosciences), and granzyme B (clone GB11; BD Biosciences) antibodies for 30 min at 4 °C. Flow cytometry data were acquired using the BD Accuri C6 and BD Accuri C6 Plus (BD Biosciences). The supernatants of these samples were subjected to Granzyme B secretion analysis using an ELISA kit (R&D Systems). Following detection of the activated cell populations in certain samples, a microfluidics-based flow cell sorting instrument (WOLF NanoCellect, San Diego, USA) was used to sort CD3⁺CD137⁺ cells from the co- culture of tumor spheres and TILs. The RNA of the sorted cells was extracted using a Nucleospin RNA XS kit (Macherey-Nagel), and its concentration was measured using QuBit. Libraries were prepared for TCR profiling using a SMARTer Human TCR a/b Profiling Kit v2 (Cat. No. 63478) and subsequent NGS library preparation and bead purification, Tapestation 4200 (Agilent) was used to confirm sample quality control with automated electrophoresis.

Statistical analysis for spheroid assays

Statistical tests for the spheroid assays were done using two-tailed unpaired t-tests, with p-values represented as * for p < 0.05 and ** for p < 0.01; all error bars represent standard error (SEM), unless otherwise stated.

Single-cell RNA-seq data analysis

Library preparation

5’ GEX and TCR sequencing: Single cell suspensions prepared from cryopreserved uveal melanoma biopsies were subjected to dead-cell removal using beads (Miltenyi Biotech) before quality and viability validation and loading on Chromium instrument (10x Genomics). Single cell 5’ gene expression and TCR libraries were prepared using vendor protocol with the Next GEM Single Cell 5’ Kit v2 (10x Genomics, catalog number 1000263), Library Construction Kit (Catalog number 1000190) Single Cell Human TCR Amplification Kit (Catalog number 1000252) and Dual Index kit TT Set A (Catalog number 1000215). Libraries were sequenced for quality check on iSeq100 (Illumina) before deep sequencing on NovaSeq 6000 (Illumina) in the required format to yield at least 20,000 reads per cell for Gene expression and 5,000 reads per cell for TCR libraries.

Alignment and gene expression quantification

Paired single-cell RNA- and TCR sequencing data from biopsies were aligned to the 10x Genomics-provided GRCh38 reference genome (v. 2020-A) and V(D)J Reference (v. 7.0.0) using the Cell Ranger pipeline (v. 7.0.1) and the “multi” command, with both “Gene Expression” and “VDJ-T” analyses activated, as well as with the setting “check-library-compatibility,false”.

Removal of ambient RNA

To remove ambient RNA, CellBender (v. 0.2.0)⁶⁵ was used on the output from Cell Ranger, with the parameters “--fpr 0.01 --epocs 150” together with expected and total cell numbers estimated from Cell Ranger quality control plots, as described in the instructions at https://cellbender.readthedocs.io/en/latest/usage/index.html (accessed 15 Sept 2022). The values for the parameters “--low-count-threshold” and “--learning-rate” were changed from their defaults for some problematic samples and iteratively adapted as needed after inspecting CellBender output quality control plots, as described in the documentation.

Separation of human and mouse reads in PDX scRNA-seq data

To separate reads of mouse and human origin in the scRNA-seq data from PDX samples, XenoCell (v. 1.0)^66,67 was used. As input, compressed fastq files, merged across lanes, were provided, together with reference genomes for human and mouse (the same genomes as for Cell Ranger). An index was created with the command “xenocell.py generate_index”, with the parameters “--threads 1 --memory 24”. Reads were classified according to host (mouse) or graft (human) with the command “xenocell.py classify_reads” and the parameters “-- barcode_start 1 --barcode_length 16 --threads 1 -memory 24 --compression_level 1”. Reads belonging to human and mouse, respectively were then extracted into new fastq files using the command “xenocell.py extract_cellular_barcodes” with the same values for barcode_start, barcode_length, threads, memory and compression_level parameters as in the previous step, but with “--lower_threshold 0 --upper_threshold 0.1” for graft and “--lower_threshold 0.9 --upper_threshold 1.0” for host. The workflow was implemented in a Nextflow script available at (https://github.com/jowkar/xenocell_nextflow). The new fastq files were then used for separate Cell Ranger runs for human and mouse reads for a given sample, as described above for human biopsy samples, except with the corresponding 10x Genomics-provided mouse reference genome instead for reads determined to be of mouse origin.

Quality control and cell filtering

Gene counts and assembled TCR sequences where imported to R with Seurat⁶⁸. An initial object was created and normalized with the NormalizeData function of the same package. Samples were integrated with fastMNN⁶⁹ (batchelor package v. 1.10.0) via the convenience function RunFastMNN from the SeuratWrappers (v. 0.3.0) package (parameters: “features = 4000, auto.merge = T, d = 150, k = 15”). Assembled and annotated TCR sequences from Cell Ranger were further added to this object using djvdj (v. 0.0.0.9000, https://github.com/rnabioco/djvdj). Cells with < 200 reads were then removed.

A first set of potential doublets were identified using scDblFinder (v. 1.8.0), using default parameters. In addition to this, cells expressing conflicting lineage-associated markers were also marked as doublets. Incompatible marker sets were defined as any among CD3D, CD4, CD8A, CD8B or with assembled TCR receptor expressed together with any of MLANA, PMEL, TYR, HBA1, HBA2 or HBB; any cell with assembled TCR receptor expressed together with any of CD14, CD19, MS4A1, or JCHAIN; or any among MLANA, PMEL or TYR expressed together with CD14, CD19, MS4A1, JCHAIN or NCR1. For the purpose of further refining doublet detection, an initial clustering was performed using the Seurat functions FindNeighbors (parameters: ‘reduction = “mnn”, dims = 1:150, k.param = 5’) and FindClusters (“resolution = 5”). For each cluster, a two-sided Fisher’s exact test was carried out, assessing whether that cluster was associated with doublets. p-values were adjusted for multiple testing with Benjamini-Hochberg correction. In addition, each cluster was also assessed for higher than expected UMI counts per cell. Any cluster with Fisher test odds ratio < 1 and q < 0.05, percent doublets greated or equal to the 95% percentile of doublet percentages, or UMI counts above the 95% percentile were further nominated as a doublet cluster to be removed. These thresholds were chosen in conjunction with visual inspection of plots regarding how these candidate clusters associated with other metrics of poor quality, such as proximity in UMAPs to cells expressing incompatible lineage markers and higher/lower than usual ribosomal or mitochondrial read percentage.

Global thresholds to keep cells fulfilling the following criteria were set: number of detected genes > 200, UMI counts > 500, ribosomal read percentage < 40% and percentage of genes in the top 100 expressed genes < 80%. However, further specific filtering was made for two cell types, erythrocytes and plasma cells, known to have generally lower RNA content. Since common global thresholds on RNA and number of detected genes may filter out these cell types entirely, lower thresholds were motivated specifically for these. To establish suitable thresholds, candidate members of these cell types were first identified based on the distinctly expressed marker genes HLA1 (erythrocytes) and JCHAIN (plasma cells, with IL3RA and LILRA4 additionally demanded to be zero), respectively. Cells in each respective group (without having removed candidate doublets) were clustered based on Spearman correlation coefficients with the pheatmap function in the R package of the same name (default parameters) and visualised together. In both cases, cells were divided into clear clusters, some of which had predicted doublets were overrepresented. The latter clusters also tended to have higher than expected RNA content. Cluster without overrepresented candidate doublets were retained and the rest discarded from the dataset. Thresholds to keep cells with UMI counts and number of detected genes greater than the first quartile were defined for each of the two cell types based on cells in clusters not associating with doublets.

Further to this, additional quality filtering thresholds for remaining cells were determined adaptively using miQC (v. 1.2.0)⁷⁰, which first creates mixture models to assess the distribution of number of genes detected versus percent mitochondrial reads in each sample and then filters cells based on a posterior probability of being compromised derived from this distribution. For a few samples that failed to converge with the default model, a one-dimensional instead of mixture model was instead used. For all samples, the parameters “posterior_cutoff = 0.9” and “keep_all_below_boundary = TRUE” were used to filter out low-quality cells.

Sample integration, unsupervised clustering and cell type annotation

After removing low-quality cells and doublets, all samples were reintegrated using fastMNN, with the same parameters as described above. A UMAP embedding was created using the Seurat function RunUMAP (parameters: reduction = “mnn”, dims = 1:150). Unsupervised clustering was performed with FindNeighbors (reduction = “mnn”, dims=1:150) and FindClusters (resolution = 2). Genes differing between clusters were assessed using the FindAllMarkers function (default parameters). Clusters were annotated based on a combination of generally known cell type markers (CD3D and CD3G for T cells, CD8A for CD8⁺ T cells, CD4 for CD4⁺ T cells, FOXP3 for regulatory T cells, NCAM1 and NCR1 for NK cells, CD19 and MS4A1 for B cells, HBA1, HBA2 and HBB for erythrocytes, CD14 for general monocytic cells, PMEL, MLANA, and TYR for melanocytic cells, as well as MKI67 for actively cycling cells) and markers compiled from literature^28,71–74. CD8⁺ T cells where determined and studied in greater detail to identify further subgroups. After subsetting the main Seurat object to only contain CD8⁺ T cells, these were subclustered by rerunning FindNeighbors (‘reduction = “mnn”, k.param = 15, dims=1:150’) and FindClusters (“algorithm = 1, resolution = 1.3”). A new UMAP was created using RunUMAP (‘reduction = “mnn”, dims=1:150, min.dist = 0.05, n.neighbors = 15’). 16 clusters were obtained, among which some were merged based on similar expression of key marker genes. The new CD8⁺ T cell clusters were annotated based on a combination of literature-derived markers^28,75–77 for CD8⁺ T cell phenotypes as well as genes found to distinguish the clusters that could be related to specific phenotypes.

Processing of scRNA-seq from TIL cultures

Single-cell RNA-seq data of TIL cultures from a previous study⁴¹ were obtained, loaded into a Seurat object in R and normalized with the NormalizeData function (default parameters). Variable features were determined using FindVariableFeatures (selection.method = “vst”, nfeatures = 4000). Cell cycle scores were then calculated using CellCycleScoring (default provided gene sets) and data scaled using ScaleData (vars.to.regress = c(“S.Score”, “G2M.Score”)). Prinicipal components were calculated using the RunPCA function (npcs = 30) and a UMAP embedding calculated using the RunUMAP function (dims = 1:30). These steps were performed on each sample separately, without sample integration.

GLIPH2 analysis

TCR sequences were clustered into inferred specificity groups using GLIPH2⁴³. As input, unique TCR sequences and their counts from annotated CD8⁺ T cells were given, together with previously published HLA genotypes from bulk RNA and WGS⁴¹. HLA genotypes for two additional samples were inferred from previously published exome data⁷⁸ using nf-core/hlatyping (v2.0.0)⁷⁹, with the parameter “--genome GRCh37”. Clonotypes with more than one TCR alpha chain were excluded from analysis, as these are not supported by GLIPH2. Parameters provided were the option to use the v2 human reference for CD8 T cells and with “all_aa_interchangeable” set to Yes. The tool was run from the web interface and according to instructions at http://50.255.35.37:8080/ (accessed 21 Dec 2022). High-confidence TCR clusters were considered those with V-gene bias score < 0.05, at least 3 unique TCRs and at least 3 unique individuals contributing TCRs.

Copy number analysis

To determine copy number changes in UM cells from scRNA-seq data, the inferCNV R package (v. 1.10.1)⁸⁰ was used according to instructions provided at https://github.com/broadinstitute/infercnv/wiki (accessed 2 Nov 2022). A genomic position file was created from the reference annotation provided with Cell Ranger (v. 2020-A) using the inferCNV helper script gtf_to_position_file.py. All cells annotated as UM were marked as “malignant” and remaining cells were used as a reference diploid background as input to inferCNV, together with the following parameters: ’cutoff=0.1, cluster_by_groups=TRUE, denoise=TRUE, HMM=FALSE, num_threads = 1, analysis_mode=“samples”, output_format = “pdf”’.

Differential expression analysis

Differentially expressed genes between experimental conditions in scRNA-seq data from PDX samples were assessed using a pseudo-bulk approach together. First, the all human CD8+ T cells for each sample were used to create gene subset expression matrix. All cells in a given sample were summed on a per-gene basis to create pseudo-bulk samples, imitating a bulk RNA-seq experiment. These samples were normalized and analysed for differential expression between the conditions represented by the samples using DESeq2⁸¹, with the parameters “test = “LRT”, reduced = ∼ 1”, relative to a full model specification encoding the two experimental conditions to be compared. For each pairwise comparison of conditions, only samples in those two conditions were included in the DESeq2 main object. To extract a table of results, the “results” function from the same R package was used, with the parameter “alpha = 0.05”. Significant genes were visualized with the pheatmap function from the R package of the same name (v. 1.0.12), with the parameter ’scale = “row”’.

scRNA and TCR-seq analysis for UM1 sample

FASTQs files lane-1 and lane-2 were concatenated and used as an input for the XenoCell (v.1.0) tool to split scRNA-seq data on host (mm) and graft (hs). XenoCell classify_reads function was used to calculate the percentage of graft- and host-specific reads for each cellular barcode. Then, extract_cellular_barcodes function was used to extract cellular barcodes for the host (mm), and the same function was used to extract cellular barcodes for the graft (hs). The feature-barcode matrix for graft (hs) was generated using cellranger multi (v.7.0.1) pipeline. Parameters --localcores and --localmem were set to 32 and 240 respectfully and human reference genome GRCh38 (v.2020-A) and V(D)J reference (v.7.0.0) were used. Next, cellbender remove- background CellBender (v.0.2.0) function was used to remove ambient RNA. For each sample cellbender parameters were adjusted as follows --expected-cells (25L: 7573; 25S: 2274; 39L: 7973; 39S: 10527; 45L: 8438; 45S: 10447; 55L: 6318; 55S: 4143), --total-droplets-included (25L: 17000; 25S: 10000; 39L: 12000; 39S: 17000; 45L: 16000; 45S: 19000; 55L: 15000; 55S: 13000), --fpr (25L: 0.01; 25S: 0.01; 39L: 0.01; 39S: 0.01; 45L: 0.01; 45S: 0.01; 55L: 0.01; 55S: 0.01), --epochs (25L: 150; 25S: 150; 39L: 200; 39S: 150; 45L: 150; 45S: 200; 55L: 150; 55S: 150), --low-count-threshold (25L: 3; 25S: 5; 39L: 3; 39S: 3; 45L: 3; 45S: 5; 55L: 5; 55S: 5), --learning-rate (25L: 0.0000125; 25S: 0.00005; 39L: 0.00000625; 39S: 0.0000125; 45L: 0.00000625; 45S: 0.00000625; 55L: 0.00005; 55S: 0.00005).

R (v.4.3.2) was used for downstream analysis. Data generated by CellBender was parsed using Read_CellBender_h5_Multi_Directory function from scCustomize (v.2.0.1) package. Seurat object was created using CreateSeuratObject function from Seurat (v. 5.0.1) package with parameters min.cells and min.features set to 3 and 200 respectfully. This resulted in the Seurat object with 23,884 genes across 64,091 cells. V(D)J data was added to the Seurat object using import_vdj function from djvdj (v. 0.1.0) library. Doublets were removed using scDblFinder function from scDblFinder (v. 1.16.0) package that resulted in 57,200 cells. Cells with a low coverage were removed before using scDblFinder (nCount <200, https://bioconductor.org/packages/release/bioc/vignettes/scDblFinder/inst/doc/scDblFinder.html). QCs filtering was performed to remove dying/dead, any remaining empty droplets or multiples cells. Cells with a high gene deletion number could indicate doublets or multiplets, thus, cells that have unique gene counts over 6,500 or less than 400 were removed. Cells with a total read count (or UMI) less than 500 and more than 25,000 were filtered out. Finally, cells with mitochondrial counts more than 20% and ribosomal counts more than 15% were removed. Next, the data was normalized, scaled, and dimensionality reduction applied using NormalizeData, FindVariableFeatures (nfeatures = 2000), ScaleData, RunPCA. This step followed standardized Seurat workflow (“Data analysis workflow Integrative analysis in Seurat v5”, https://satijalab.org/seurat/articles/seurat5_integration, accessed 2024). Integration was performed using IntegrateLayers function from Seurat (v. 5.0.1) package with parameters set to method = CCAIntegration, orig.reduction = “pca”, new.reduction = “integrated.cca”. FindNeighbors (reduction = “integrated.cca”, dims = 1:70) and FindClusters (resolution =c(0.1 to 2.0) , cluster.name = “cca_clusters”) functions from Seurat (v. 5.0.1) package were used and resolution 0.6 (15 clusters) was selected. FindAllMarkers (logfc.threshold = 0.3, min.pct = 0.3, only.pos = TRUE, test.use = “wilcox”) function from Seurat (v. 5.0.1) package was used to find differentially expressed genes. To visualize results DimPlot, VlnPlot, FeaturePlot and DotPlot Seurat (v. 5.0.1) functions were used. For plotting the volcano plot we used EnhancedVolcano (v.1.20.0; accessible from https://github.com/kevinblighe/EnhancedVolcano), for plotting a bar plot ggplot2 (v.3.5.0) was used. To pseudobulk, AggregateExpression function from Seurat (v. 5.0.1) package was used. The gene list was compiled using marker genes coming from differential expression analysis and from the literature. Results were visualized as a heatmap using plotGroupedHeatmap function from the scater (1.30.1) package. To find clonotype abundance calc_frequency function from djvdj (v. 0.1.0) package was used, to visualize results plot_clone_frequency function was used from the same package.

TCR sequence analysis

Bulk TCR sequencing data analysis

TCR sequences sequenced with the Takara SMARTer Human TCR a/b Profiling Kit v2 were aligned and assembled using the open source Cogent NGS Immune Profiler software (v1.0) provided by the same company, using the parameters “-r TCR” and “-t Both”.

Matching of TCR receptors between samples

TCRs from single-cell RNA-seq samples and/or bulk TCR-seq samples were matched based on identical combinations of CDR3β, Vb and Jb. This was the minimal common combination of criteria available on which to consistently compare TCRs across all TCR-sequenced samples, since bulk the TCR-seq protocol does not provide information in pairing of alpha and beta chains.

Matching of TCR receptors to public databases

To match TCRs to public databases with information on possible antigens they might recognize, TCRMatch (v. 1.0.2)⁴² was used (parameter: -s 0.97) together with a modified version of the reference database from IEDB⁸² provided with the tool, which we expanded with additional information from VDJdb⁸³, McPAS-TCR⁸⁴ and TCR3d⁸⁵ (using the datasets for both cancer and virus-targeting TCRs). In addition, this database file was further modified to replace internal commas in antigen names, since the default formatting conflicts with the use of comma as the main delimiter in the database file. This facilitated easier analysis of the results. TCR matches were visualized with Sankey diagrams using the ggsankey R package (v. 0.0.99999, https://github.com/davidsjoberg/ggsankey).

Overrepresentation analysis for TCRs in CD8⁺ T cell clusters

To assess whether the 41BB- or MART1-reactivity experiment would enrich for clones mapping to any phenotype relative to cultured TILs in general, a binomial test was used. The test considers the observed fraction of cells with matches to the reactivity screen versus an expected probability. Given that cells from the reactivity experiment were taken from a TIL culture, we wanted to assess whether these had a meaningfully different distribution compared to cultured TILs in general in regards to cluster representation. Therefore, the expected probability was estimated as the fraction of TILs that match any given cluster. Thus, these two frequencies were compared with one two-sided binomial test (binom.test in R) per cluster, and p-values further adjusted for multiple testing with Bonferroni correction.

Data access statement

The single-cell sequencing data has been uploaded to European Genome-Phenome Archive and is available under controlled access after publication. All remaining data are within the paper and its Supporting Information files or available on request.

Author contributions

JK: conducted all bioinformatic analysis except for the UM1 liver metastasis model data, generated figures and wrote several sections of the paper; VS: conducted experiments, generated figures and wrote some sections of the paper; ROB: coordinated patient consents, surgeries and collection of liver metastases of uveal melanoma; MI: conducted single-cell sequencing; SA: performed cell sorting and flow analyses; SS: generated PDX models and developed new transplantation and imaging techniques; IK: conducted bioinformatic analysis for the liver metastatic model of the UM1 sample; AS: managed single-cell sequencing workflow; LN: supervised the study, coordinated patient consents and sample collection; LMN: supervised, managed the biobank and sequencing, generated TILs, managed, conceived and visualized data from animal experiments; JAN: conceived and supervised the study, wrote the paper and generated figures; all authors, read, edited and approved the manuscript.

Acknowledgements

We thank Carina Karlsson for providing technical support. Grant support came from Cancerfonden (to J.A.N. and R.O.B), Familjen Erling Persson (to J.A.N.), Knut and Alice Wallenberg Foundation (to J.A.N. and R.O.B), Vetenskapsrådet (to J.A.N. and R.O.B.), Sjöbergstiftelsen (to J.A.N.), BioCARE Strategic grants (to J.A.N.), Lion’s Cancerfond Väst (to S.A. and L.N), Västra Götaland Regionen ALF grant (to J.A.N., R.O.B. and L.N.), Assar Gabrielsson fond (to V.S.), Gustaf V Jubileumsklinikens forskningsfond (to L.N.) and Wilhelm & Martina Lundgrens Vetenskapsfond (to J.K.). J.A.N. is the Inaugural Chair of Melanoma Discovery, primarily supported by donations to Perkins from family and friends of Scott Kirkbride, the MACA Cancer 200 ride for cancer research, the HBF Foundation and Perpetual. The Genomics WA facility (A.S.) is supported by BioPlatforms Australia, State Government Western Australia, Australian Cancer Research Foundation, Cancer Research Trust, Harry Perkins Institute of Medical Research, Telethon Kids Institute and the University of Western Australia. We also thank the Pawsey Supercomputing Centre for providing computational resources for this project.

Supplementary Figure and Table Legends

Supplementary Figure 1.

Supplementary Figure 2.

Supplementary Figure 3.

Supplementary Figure 4.

Supplementary Figure 5.

Supplementary Figure 6.

Supplementary Figure 7.

Supplementary Table 1. Statistically identified marker genes for each CD8⁺ T cell cluster in Fig. 2a., using the FindMarkers function in the Seurat R package.

Supplementary Table 2. Complete statistics for the analysis of cluster representation among reactivity-screened TCRs matched to biopsy CD8⁺ T cells, as compared to representation among TILs matched to the same clusters (Fig. 4d).

Supplementary Table 3. Complete statistical output from the GLIPH2 analysis referred to in Fig. 4h-i and Supplementary Fig. 4i.

Supplementary Table 4. Complete statistics for the differential expression analyses in Fig. 5k and Supplementary Fig. 5j.

Supplementary Table 5. Differentially expressed genes in different clusters of the UMAP in Fig. 5.

Supplementary Table 6. References for marker gene selection and annotation of clusters in in Fig. 5.