Integrative analysis of neuroblastoma by single-cell RNA sequencing identifies the NECTIN2-TIGIT axis as a target for immunotherapy

Summary Pediatric patients with high-risk neuroblastoma have poor survival rates and urgently need more effective treatment options with less side effects. Since novel and improved immunotherapies may fill this need, we dissect the immunoregulatory interactions in neuroblastoma by single-cell RNA-sequencing of 24 tumors (10 pre- and 14 post-chemotherapy, including 5 pairs) to identify strategies for optimizing immunotherapy efficacy. Neuroblastomas are infiltrated by natural killer (NK), T and B cells, and immunosuppressive myeloid populations. NK cells show reduced cytotoxicity and T cells have a dysfunctional profile. Interaction analysis reveals a vast immunoregulatory network and identifies NECTIN2-TIGIT as a crucial immune checkpoint. Combined blockade of TIGIT and PD-L1 significantly reduces neuroblastoma growth, with complete responses (CR) in vivo. Moreover, addition of TIGIT+PD-L1 blockade to standard relapse treatment in a chemotherapy-resistant Th-ALKF1174L/MYCN 129/SvJ syngeneic model induces CR. In conclusion, our integrative analysis provides promising targets and a rationale for immunotherapeutic combination strategies.


INTRODUCTION
Immunotherapy has revolutionized cancer treatment in adults and holds great promise for pediatric solid tumors. 1 Its potential is exemplified by the increased survival of patients with high-risk neuroblastoma following implementation of anti-GD2 antibody therapy into standard care. 2,3Neuroblastoma, the most common extracranial pediatric solid tumor, accounts for 10% of pediatric cancer-related deaths. 4Patients are stratified into risk groups based on disease presentation and genomic alterations like MYCN amplification, which is an important driver of poor prognosis. 5High-risk neuroblastoma patients receive an intense multimodal treatment regimen, consisting of induction chemotherapy, surgical tumor resection, consolidation with high-dose chemotherapy followed by autologous stem cell transplantation, radiotherapy and anti-GD2 immunotherapy.Introduction of anti-GD2 has improved high-risk neuroblastoma event-free survival rates by $15%. 2,3Still, overall 5-year survival rates are below 60%, particularly due to the high relapse rate. 3Chemoimmunotherapy with temozolomide, irinotecan and anti-GD2 is likely to be selected as the backbone treatment for relapse/refractory neuroblastoma patients in Europe. 6Due to the intense treatment, the majority of survivors suffer of debilitating (long-term) side effects. 7Taken together, there is an urgent need for more effective treatments with less side effects.
The success of anti-GD2 therapy has provided a clear rationale for immunotherapy in neuroblastoma treatment.Yet, as the survival benefit is still modest, improving immunotherapy efficacy will be crucial.T cells and natural killer (NK) cells are considered essential effectors in the context of immunotherapy, but their function is compromised in many cancers. 8,9Efforts to improve immunotherapy efficacy focus on reinvigorating T and NK cell function, e.g., with immune checkpoint inhibition or by introducing cellular immunotherapies like chimeric antigen receptor (CAR) T cells. 1,102][13][14][15][16][17][18][19][20] The overall limited success may be due to intra-tumoral factors hampering efficacy, such as immunosuppressive cells and a plethora of immune checkpoints. 10,21Efforts to optimize immunotherapies for neuroblastoma are currently hampered by a lack of information on these intra-tumoral factors.The detailed composition and function of T and NK cells in the neuroblastoma tumor-microenvironment (TME) is still largely unexplored and only few immunosuppressive factors have so far been identified. 22To guide further immunotherapy development-and enhancement-it is essential to gain better insights into neuroblastoma's immune environment and the immunoregulatory mechanisms at play.Moreover, such insights may fuel an educated approach to design combination immunotherapy strategies overcoming immune resistance. 22o unravel the immune environment of neuroblastoma and identify targets for immunotherapy enhancement, we generated a single-cell transcriptomic atlas of 24 tumors, with samples taken before and after induction chemotherapy.We provide a detailed view of neuroblastoma's immune landscape pre-and post-chemotherapy and reveal a vast immunoregulatory interaction network.Subsequent functional validation experiments identified the NECTIN2-TIGIT axis as promising target for neuroblastoma immunotherapy.

RESULTS
The single-cell landscape of neuroblastoma To provide an in-depth view of neuroblastoma's TME, tumor samples were analyzed by single-cell RNA sequencing.Fresh material was processed from 24 tumors of 19 patients, 10 taken pre-treatment and 14 after induction chemotherapy, including five paired samples.Seventeen patients had high-risk neuroblastoma and two had intermediate risk neuroblastoma.Six out of 19 patients had MYCN amplified (MYCN-A) tumors (Figure 1A; Table S1).Single-cell RNA sequencing using the Cel-Seq2 protocol yielded 22,418 high-quality cells.We identified 24 cell clusters, which were annotated as five main cell types, i.e., tumor, immune cells, endothelium, mesenchyme, and Schwann cell precursor cells (Figures 1B, 1C, and S1A), confirming previous findings. 23,24The three tumor clusters highly expressed CDK4, CCND1, and MYCN, respectively (Figure 1C).Cluster 1 contained CDK4 amplified cells and cluster 3 contained MYCN-A cells (Table S1).In addition, the tumor clusters expressed previously described neuroblastoma markers genes, including CHGB, PHOX2B, and TH (Figures 1D and S1B). 25Tumor cell identity was confirmed by copy number variation inference, showing typical chromosomal aberrations (e.g., 1p loss, 11q loss, and 17q gain; Figure S1C). 25LA class I expression was significantly lower in tumor than nontumor cells, particularly in those with high MYCN expression (Figures S1D-S1F), highlighting neuroblastoma's low immunogenicity. 26Nonetheless, the cellular (immune) fractions of MYCN-A tumors did not significantly differ from MYCN-non amplified tumors (Figure S1G).Distinct changes in the cellular composition were observed upon induction chemotherapy (Figures 1E-1G).The tumor content decreased from 47% (range 9-96%) pre-treatment to 3% (0-20%) post-treatment, while the mesenchymal fraction increased from 11% (0-43%) to 49% (1-87%) (both p < 0.05, Figures 1F and 1G).The mesenchymal fraction consisted of (cancer-associated) fibroblasts (CAF) with phenotypes ranging from myofibroblasts to inflammatory/adipogenic fibroblasts, one population enriched for follicular dendritic cells (FDC) derived from lymph node biopsies (C8), and one population of healthy adrenal cells (C9; Figures S1H-S1J).CAFs displayed functional and compositional differences before and after treatment (Fig- ures S1I-S1K).C0 and C8, containing myofibroblasts and inflammatory CAF/FDC, were the most prominent CAF clusters pretreatment, while clusters C2, C5, and C6 were more enriched post-treatment.The other main cell fractions, including immune cells, did not change upon treatment (Figure 1F).

Neuroblastoma is populated by immunosuppressive macrophages
To unravel the composition and function of neuroblastoma-infiltrating immune cells we performed an in-depth analysis of the PTPRC + (CD45 + ) immune clusters.Among a total of 6,012 immune cells, we identified myeloid and lymphoid populations (Figures 2A,  S2A, and S2B).The myeloid cells, key regulators of anti-tumor immunity due to their roles in antigen presentation, T cell polarization and immunosuppression, 27 consisted of mast cells, plasmacytoid dendritic cells (pDC), conventional DC (cDC), S100A8/A9 hi undifferentiated monocytes (Mo), and four differentiated macrophage populations (IL10 hi , CCL2 hi , APO hi , and MAF hi macrophages [M4; Figure 2B]).We confirmed their presence by flow cytometry and established their similarity with previously identified populations in neuroblastoma (Figures S2C-S2E). 23,24DC, and IL10 hi , APO hi and MAF hi M4 had the highest expression of HLA and costimulatory molecules (p < 0.0001 versus S100 hi Mo), rendering them the most potent interaction partners for T cell co-stimulation (Figures 2C and S2F).Compared to the undifferentiated S100 hi Mo, all macrophage populations displayed an M2-like signature, associated with immunosuppressive and pro-tumorigenic properties (p < 0.0001; Figure 2D), 28,29 which was confirmed by protein expression of M2-associated marker CD163 (Figure 2E).The IL10 hi and CCL2 hi M4 populations had a mixed profile, expressing M1-like (pro-inflammatory) features next to M2-like characteristics (p < 0.0001; Figure 2D).1][32] These results indicate, and confirm previous reports, 23,24 that macrophages in neuroblastoma have immunosuppressive features which may tune lymphoid responses.

Lymphocyte subsets in neuroblastoma display features of dysfunctionality
The lymphoid compartment consisted of NK cells, T cells, and (plasma) B cells (Figure 3A).NK cells expressed genes encoding cytotoxic molecules (GNLY, GZMA, GZMB, and PRF1) and gdT cells expressed their T cell receptors (TCR; TRDC, TRGC1/2).Among abT cells, we identified CD8 + , CD4 + , CD4 + FOXP3 + regulatory T cells (Tregs), and naive(-like) T cells.These populations were largely equivalent to previously identified clusters in neuroblastoma (Figure S3A). 23,24Histological assessment revealed that T cells were present in tumor-rich areas (Figure S3B).CD8 + T cells expressed CD8A, cytotoxic effector molecules  (PRF1, GZMA, and GZMB), and significantly increased levels of PDCD1 (encoding PD-1) and LAG3 compared to other T/NK cell clusters (Figures 3A and 3B).Since PD-1 and LAG-3 are typically regarded as markers of T cell dysfunction/exhaustion, their expression indicated the presence of a dysfunctional cell fraction among CD8 + cells. 33,34Tregs expressed high levels of their signature genes FOXP3, CTLA4, TIGIT, and IL2RA, in addition to transcription factors IKZF2, BATF, PRDM1, and MAF and checkpoint receptors TNFRSF1B (TNFR2), TNFRSF4 (OX-40), and TNFRSF18 (GITR) (Figures 3A, 3C, and S3C).][37][38][39] Among non-Treg CD4 + T cells we identified two distinct subsets with functional differences.One subset highly expressed DUSP4, a signal repressor preventing T cells from over-activation, which has been related to premature T cell aging causing senescence/ exhaustion (Figure 3D). 40DUSP4 hi CD4 indeed had a mixed profile of activation and regulation, with on the one hand TCR signaling and expression of proliferation marker MKI67, and on the other hand expression of co-inhibitory receptors PDCD1, TIGIT, CTLA4 and high PD-1 signaling (Figures 3D and S3D-S3F).Expression of these co-inhibitory receptors by non-Treg T cells is associated with T cell dysfunction. 34,41Indeed, the DUSP4 hi CD4 + cluster had significantly increased expression of a published signature for dysfunctional tumor-infiltrating T cells (Figure 3E). 33In contrast, TPT1 hi CD4 + T cells had higher translational activity and expressed receptors IL7R and CCR7, cytotoxic markers (KLRB1, GNLY), as well as CCR6 and RORA, associated with Th17 polarization (Figures 3D and S3F).We confirmed the presence of these distinct CD4 + T cells in neuroblastoma by flow cytometry (Figure 3F), including high expression of co-inhibitory receptors and Ki-67 in the DUSP4 hi cluster (Figures S3G and S3H).These results suggest that DUSP4 hi CD4 contained a higher proportion of tumorreactive cells, whereas TPT1 hi CD4 likely consisted of tissue-resident bystander cells (Figure S3I).This premise was further supported by gene set enrichment analysis (GSEA) utilizing published gene signatures specific for either tumor-infiltrating T cells or blood/normal tissue-associated T cells. 42GSEA showed that DUSP4 hi CD4 corresponded largely with previously identified tumor-infiltrating T cell clusters-including those with exhausted profiles-whereas TPT1 hi CD4 corresponded more with blood/ normal tissue-associated Th17 cell clusters (Figure 3G).Bystander TPT1 hi CD4 may have been attracted to tumors by CCR6 ligand CCL20, which was highly expressed by IL10 hi M4 (Figure S2F).Taken together, the neuroblastoma immune environment is characterized by lymphocyte subsets with various features of dysfunctionality, and highly immunosuppressive effector Tregs.
NK cells in pre-treatment neuroblastoma tumors are dysfunctional Nearly all annotated immune subsets were detected in each individual tumor sample, albeit in varying proportions (Figures S4A  and S4B).To assess the effect of induction chemotherapy on tumor immunity we compared the immune composition before (at diagnosis) and after induction chemotherapy (at surgical resection) (Figure 4A).While the total proportion of immune cells was  constant (Figure 1F), we observed a trend of a decreased lymphoid/myeloid ratio upon treatment, suggesting a reduction of lymphoid or an influx/expansion of myeloid cells (p = 0.13; Figures 4B and S4A-S4C).The strongly M2-differentiated CCL2 hi M4 significantly increased after treatment, while B cells decreased (both p < 0.05; Figures 4C, 4D, S4D, and S4E).We also observed a trend of NK cell reduction after treatment (p < 0.1; Figures 4D and S4E).Since NK cells are considered essential cytotoxic effectors in neuroblastoma, 43 we assessed whether their functionality changed upon treatment, and compared them to reference NK cells from healthy donor blood (Figures S4F and S4G).Since NK (G) GSEA comparing neuroblastoma immune cell clusters (in rows) to previously published gene signatures which identify T cell clusters isolated from either blood, healthy tissues (''normal'') or tumors (in columns). 42The color scale indicates for each neuroblastoma cluster the degree of similarity (NES) with published signatures.Also see Figure S3.(legend continued on next page) cell function could be affected by micro-environmental effects induced by tumor cells, which were most abundant pre-treatment (Figure 1F), we primarily focused on NK cell function pre-treatment.
Compared to NK cells post-treatment or reference NK cells, NK cells in pre-treatment tumors had significantly reduced expression of essential cytotoxic effector genes (GNLY, GZMB, PRF1; Figures 4E, 4F, and S4H).We confirmed low granzyme B and perforin expression, suggesting an impaired cytotoxic function, by flow cytometry (Figure 4G).GSEA further confirmed their reduced cytotoxicity and indicated a dysfunctional immature, resting state in pre-treatment NK cells (Figure S4I), as also observed in other tumor types. 44Signaling by TGF-b1, a well-known immunosuppressive factor, correlated negatively with NK cell cytotoxicity (Figure 4H). 45TGF-b downstream signaling was observed in all TME cell types and correlated strongly with TGF-b receptor expression, supporting the presence of active TGF-b in neuroblastoma (Figures S4J and S4K).Intriguingly, NK cells themselves expressed the highest levels of TGF-b1 (Figure S4J).Addition of recombinant TGF-b1 to in vitro co-cultures of neuroblastoma tumoroids and primed healthy donor NK cells significantly reduced tumor killing and granzyme B production by NK cells, whereas blockade of TGF-b increased cytotoxicity (Figures 4I and 4J).These observations strongly implicate TGF-b1 signaling in NK cell dysfunction in neuroblastoma.
Since NK cell activity is additionally regulated by a delicate balance of activating and inhibitory receptors, we assessed expression of these receptors and evaluated their correlation with NK cell cytotoxicity (Figure S4L).Overall, activating receptor expression was lower in tumor than blood NK cells (Figure S4M).Expression of activating receptors CD226 (DNAM-1), FCGR3A (CD16), KLRC2 (NKG2C), KLRD1 (CD94) and NCR1 (CD335) positively correlated with cytotoxic genes, while the inhibitory KLRB1 (CD161), CD96 and TIGIT negatively correlated with cytotoxic genes (Figure 4K).The activating CD226 and inhibitory CD96 and TIGIT belong to the same checkpoint receptor family competing for ligands.TIGIT/CD226 and CD96/CD226 gene expression ratios were substantially increased in tumor-infiltrating NK cells pre-treatment (Figures 4L and S4N), indicating a disturbed balance shifted toward NK cell inhibition.We confirmed significantly increased CD96 expression, and a trend of increased TIGIT expression in neuroblastoma-infiltrating NK cells by flow cytometry (Figure S4O).These results implicate TGF-b1 and the inhibitory checkpoint receptors KLRB1, TIGIT, and CD96 in NK cell dysfunction in neuroblastoma.

T cells in post-treatment neuroblastoma tumors are dysfunctional
In contrast to NK cells, T cells did not have reduced cytotoxicity pre-treatment (Figure S5A).Pathway analysis showed that gd and CD8 + T cells had higher IL-12 signaling pre-treatment than post-treatment (Figures S5B and S5C), suggesting activation.CD8 + and DUSP4 hi CD4 + T cells additionally had increased interferon and TCR signaling pre-treatment (Figures S5C and S5D), indicating that at least a fraction of those cells may be tumor-reactive, corroborating our findings in Figure 3G.All four effector T cell subsets shared pathways related to TCR signaling in the top 50 of upregulated pathways pre-treatment (Figure S4E).GSEA confirmed increased TCR signaling in all T cell subsets pre-treatment compared to post-treatment (Figure S4F).7][48][49][50] GSEA confirmed a significant enrichment of gene signatures for dysfunction in post-treatment compared to pre-treatment ab T cells, whereas ab T cells in pre-treatment samples had a higher effector signature (Figure 4N).ab T cells posttreatment were particularly enriched for signatures associated with ''terminal exhaustion'' as opposed to ''progenitor exhaustion,'' the latter being associated with some level of retained functionality and tumor control. 51However, compared to healthy donor blood-derived T cells, both pre-and post-treatment tumor-infiltrating T cells showed a significant enrichment of dysfunction/ exhaustion signatures (Figure S5I).The co-inhibitory receptors LAG3 (LAG-3), CTLA4 (CTLA-4), PDCD1 (PD-1) and HAVCR2 (TIM-3), extensively linked to T cell exhaustion/dysfunction, 33,34 were among the core enriched upregulated genes post-treatment (Figures 4O and S5J).Lower fractions of MKI67 (Ki-67) and IL2 (IL-2) expressing cells and increased fractions of TNFRSF9 (4-1BB) expressing cells confirmed reduced (proliferative) activity and suggested prolonged antigen stimulation post-treatment (Figure 4P). 53,54Taken together, these results imply that a fraction of T cells is tumor-reactive and becomes more exhausted/dysfunctional after treatment, possibly due to prolonged antigen stimulation.The high expression of co-inhibitory receptors in combination with overall low TOX/TOX2 expression (Figure 3A) however suggests potentially retained responsiveness to immune checkpoint blockade (ICB).

Immunoregulatory interactions in neuroblastoma
Lymphocyte dysfunction in tumors may result from prolonged antigen stimulation, but also from suppressive microenvironmental cues.To unravel these cues, we analyzed the interactions between T/NK cells and other cells in the tumor using the prediction tool CellChat. 52The combined population of T/NK cells (defined in Fig- ure 2A) displayed a multitude of interactions with other cells, of which those with myeloid cells were most abundant (Figures 5A, S6A, and S6B).To specifically identify therapeutically targetable immunomodulatory interactions, we focused on interactions between T/NK cells and four immunoregulatory interactions partners, i.e., myeloid cells, tumor cells, mesenchymal cells and Tregs.Interactions with these partners included, among others, interactions involved in cellular adhesion, chemoattraction and immune regulation (Figures 5B and S6A-S6F).Among these we identified a number of cellular interactions with a known immunosuppressive role (M) Venn diagram of shared upregulated genes (padj < 0.05) in CD4 and CD8 T cells post-treatment versus pre-treatment.(N) GSEA of exhaustion and effector signatures (1 + 4: GSE84105, 2 + 5, 78 3 33 ) in CD4 and CD8 T cells pre-/post-treatment. NES = normalized enrichment score.
Immunosuppressive interactions with potential for therapeutic intervention ICB has revolutionized cancer therapy and may restore activity of exhausted/dysfunctional T cells. 1 We therefore sought to identify specifically which interactions with tumor-reactive T cells were associated with disturbed T cell function, to therapeutically target exactly those interactions.We assessed the correlation between T cell dysfunction and all predicted interactions (not selected for known immunosuppressive interactions) with an unbiased approach outlined in Figure 6A, using a dysfunction score with established markers (Figures S7A-S7C). 55We focused our analysis on CD8 + and DUSP4 CD4 + T cells since these proved to contain most tumor-reactive cells and had the highest dysfunction scores (Figure S7C).The analysis revealed multiple targets which positively correlated with T cell dysfunction (Figure 6B).Among the most frequent targets and those with the highest correlations with T cell dysfunction, we identified NECTIN2, a well-known interaction partner of the co-inhibitory receptor TIGIT (Figures 6B, S7D, and S7E).NEC-TIN2 on cDC, endothelium, mesenchyme, IL10 hi M4 and APO hi M4 was predicted to interact with TIGIT on CD8 and DUSP4 hi CD4 T cells (Figures 6C, S7F, and S7G).We confirmed expression of NECTIN2 by flow cytometry (Figure 6D).Expression of TIGIT was increased in tumor-infiltrating T cells, especially in DU-SP4 hi CD4 and CD8 T cells, compared to reference T cells from blood (Figures 6E and S7H).We confirmed expression of NEC-TIN2 and TIGIT by these cells in previously published neuroblastoma datasets (Figures S7I and S7J), and confirmed the positive correlation between NECTIN2 expression and T cell dysfunction in bulk-RNAseq data of 498 neuroblastomas (SEQC cohort GSE49710; Figure 6F). 56Taken together, the NECTIN2-TIGIT axis may regulate T cell function in neuroblastoma and represents a promising target for therapeutic intervention.

Combined TIGIT/PD-L1 blockade enhances immune responses against neuroblastoma
To investigate the therapeutic potential of TIGIT blockade in neuroblastoma, we explored its functional role in tumor killing.To consider the possible benefit of ICB combination therapy-based on the abundance of immunoregulatory interactions including NK (C) Bubbleplot of selected, predicted immunosuppressive interactions with T/NK subsets.Interactions with each specific myeloid subset were evaluated and subsequently merged, with the highest probability of each interaction pair depicted in the plot.Also see Figure S6.
CD274-PDCD1 (Figure 5C) -we also tested combined TIGIT/PD-L1 blockade.PD-L1 was highly expressed by myeloid cells, particularly IL10 hi M4 (Figures S8A and S8B).Moreover, neuroblastoma-infiltrating T cells, particularly DUSP4 + CD4 T cells and CD8 T cells, displayed active downstream TIGIT and PD-1 signaling (Figure S8C).TIGIT blockade in vitro resulted in significantly enhanced killing of neuroblastoma tumoroids by immune cells, when combined with PD-L1 blockade (Figures 7A, 7B, and S8D-S8F).The rationale for this ICB combination was further substantiated by upregulation of nectin-2 and PD-L1 by tumor cells during co-culture (Figure 7C).To elucidate which immune cells contributed to the ICB-enhanced tumor-reactivity, we performed killing assays with isolated immune subsets: CD4, CD8, gdT, NK cells, and full peripheral blood mononuclear cells (PBMC) lacking each of these populations (Figure S8G).A trend of increased killing by the ICB combination was conserved in 3/3 donors in cultures with each of the four separate populations and with PBMC lacking CD4 or gdT cells.However, in cultures with PBMC lacking CD8 or NK cells, the enhanced killing efficacy was lost in 2/3 donors.CD8 T cells and NK cells might thus contribute to the ICB-induced response, and cultures/tumors without CD8 T cells and/or NK cells might benefit less from combined anti-TIGIT/PD-L1 therapy.
In vivo, combined TIGIT/PD-L1 blockade significantly improved survival and led to sustained CR in 2 out of 6 animals in two immunocompetent syngeneic murine models (N1E-115 and Neuro2a), a partial response in Neuro2a, and 1 transient CR and 1 sustained CR in a third model (N18) (Figures 7D-7G, S8H, and S8I).Moreover, the combination treatment significantly improved survival compared to anti-PD-L1 monotherapy in 2/3 models, suggesting a significant benefit of TIGIT blockade.None of the treatments led to changes in body weight or other toxicities (Figure S8J).
To understand the immunobiology of ICB in vivo we performed high-dimensional flow cytometry on the treated tumors.The percentage of CD45 + tumor-infiltrating immune cells was not affected by ICB treatment (Figure S8K).CD11b + myeloid cells were highly abundant across models and treatments (Figures 7H and S8L-S8N).Intriguingly, infiltration of CD8 + T cells, NK cells and DC correlated with a lower tumor volume in 2 out of 3 models, while infiltration with potentially (E) Flow cytometric validation of TIGIT protein expression on T cell populations infiltrating neuroblastoma, compared to reference T cells from blood (PB).TPT1 hi CD4 were gated as IL-7R hi PD-1 lo CD4 + cells and DUSP4 hi CD4 were gated as IL-7R lo PD-1 hi CD4 + cells, as shown in Figure 3F.Two-way ANOVA with Sidak's posthoc test.Mean ± SD.
(F) Correlation of NECTIN2 gene expression with dysfunction score in bulk-RNAseq dataset of SEQC cohort consisting of 498 neuroblastomas (r2.amc.nl;Tumor Neuroblastoma-SEQC-498-RPM-seqcnb1; GSE49710).Also see Figure S7.(legend continued on next page) immunosuppressive CD11b + myeloid cells and, specifically, Ly6G + myeloid-derived suppressor cells, correlated with a higher tumor volume (Figure 7I).The combination treatment significantly increased the percentage of T cells, specifically CD8 + T cells, leading to an increased CD8/Treg ratio (Figure 7J).A trend of increased NK cells and DC was also observed.On the other hand, the percentage of total CD11b + myeloid cells was decreased, with a trend of reduced Ly6G + MDSC (Figure 7J).Taken together, these results imply that TIGIT+PD-L1 blockade may engage CD8 T cells and NK cells, and may significantly modulate the immune microenvironment of neuroblastoma into a less immunosuppressive milieu.In conclusion, we identified combined TIGIT+PD-L1 blockade as a relevant intervention for neuroblastoma with therapeutic potential in vivo.
TIGIT blockade improves survival in a chemotherapyresistant neuroblastoma model Tumor relapses are the major cause of death in patients with neuroblastoma, and more effective treatments are especially urgent for these patients.We therefore investigated the added effect of TIGIT blockade in an immunologically cold, chemotherapy-resistant syngeneic model mirroring relapsed/refractory patients.The model was generated by allograft of neurospheres derived from Th-ALK F1174L /MYCN 129/SvJ transgenic spontaneous tumors subjected to repeat chemotherapy (Figures S9A-S9F). 57TIGIT ± PD-L1 blockade was added to temozolomide and irinotecan (TEM/IRI) + anti-GD2 treatment, which is favored as the standard backbone treatment for relapsed/refractory patients in Europe (Figures 8A and 8B).The combination of TEM/IRI+anti-GD2 with ICB significantly improved survival compared to the vehicle control (Figure 8C).Remarkably, TEM/IRI+anti-GD2 was already highly efficacious, and addition of ICB did not significantly improve survival compared to this backbone.Addition of TIGIT ± PD-L1 blockade did however result in three sustained (>150 days) CR (Figures 8C  and 8D).Notably, addition of anti-TIGIT to anti-PD-L1 significantly reduced tumor growth and effectuated a trend of increased survival, while PD-L1 blockade by itself did not effectuate a clear survival benefit, possibly due to relatively low expression of PD-L1 in the tumor (Figure S9G).This highlights the relevance of TIGIT as immune checkpoint in these tumors.None of the treatments led to toxicities (Figure S9H).Since the tumors in this study were relatively small (<50 mm 3 at start of treatment), we also tested TIGIT blockade in animals with larger tumors (50-200 mm 3 ; Figure S9I).In these animals, TEM/ IRI+anti-GD2 treatment alone was less efficacious, and addition of TIGIT blockade significantly improved animal survival compared to the backbone treatment alone (p = 0.0299; Figures 8E-8G, S9J, and S9K).To understand the immunobiology of the different responses in small and large tumors, we compared their tumor-microenvironments.Smaller tumors had significantly higher CD45 + immune cell infiltration (Figure 8H), resulting in significantly higher CD8 + and CD8 À CD4 À (possibly gd) T cell fractions, higher CD11b + myeloid fractions, including F4/ 80 + M4, and a trend of higher B cells and M2-like M4 as fraction of total cells in the tumor (Figures 8I, 8J, and S9L).The higher presence of effector cells might explain why the ICB treatment induced CR in a number of animals with small tumors, but not in large tumors.
Lastly, to investigate whether TIGIT ± PD-L1 blockade could add value in the context of adoptive cell therapies, we analyzed single-cell RNA sequencing data from a recent clinical trial with anti-GD2 CAR-NKT cells in neuroblastoma.This publication included analyses of isolated CAR-NKT cells pre-infusion and post-infusion. 58TIGIT expression significantly increased on CAR-NKT cells post-infusion and was highest in patients not responding to the treatment, while PDCD1 expression did not differ (Figure 8K).This finding could provide a rationale to combine adoptive cell therapies with TIGIT blockade in neuroblastoma.Taken together, TIGIT blockade may add a survival benefit in the context of chemotherapy resistance in patients with neuroblastoma, when combined with the currently favored relapse treatment protocol.

DISCUSSION
In this study we have generated a comprehensive single-cell atlas of the neuroblastoma immune environment, revealing the detailed composition and functional profile of immune cells in neuroblastoma.We exposed a vast immunoregulatory network affecting T and NK cell function and identified TGF-b1, CD161, PD-1, and TIGIT as potential targets for immunotherapy.With in vitro and in vivo studies, we demonstrated a significant benefit of anti-TIGIT/PD-L1 combination therapy in three models, a significant contribution of anti-TIGIT to this response, and a significant benefit of adding anti-TIGIT to the treatment of large chemotherapy-resistant tumors.These results nominate TIGIT(+PD-L1) blockade as a treatment option for patients with neuroblastoma, even in the context of chemotherapy resistance.
Our analyses revealed that effector lymphocytes in neuroblastoma are dysfunctional, as also observed in other tumors. 8,9,42K cells in neuroblastoma had reduced cytotoxicity, particularly pre-chemotherapy.This resting, immature profile of NK cells recapitulates their previously described state of ''arrested  (legend continued on next page) development'' in other tumors. 44NK cell dysfunction in neuroblastoma is not limited to the tumor but has also been observed in peripheral blood of patients, throughout the treatment course. 59A recently published phase II trial adding anti-GD2, considered to rely at least partly on activity of NK cells, 43 to induction chemotherapy, demonstrated improved early responses and event-free survival in comparison with historic reference cohorts. 60This added value of anti-GD2 during the induction phase-despite the dysfunctionality of NK cells in pre-treatment tumors-suggests that immunotherapy may be able to reverse or at least ameliorate the arrested state of NK cells, 60 as also previously observed in peripheral blood of patients. 59Another explanation could be that not NK cells, but rather myeloid populations are the main effector cells mediating anti-GD2 efficacy. 61hile T cells in neuroblastoma showed signs of tumor reactivity pre-treatment, they exhibited features of dysfunction post-treatment. 8,62T cell dysfunction may result from chronic stimulation of T cells and/or immunomodulatory signals. 1,63onsidering the high expression of antigen-experienced T cell marker 4-1BB post-treatment, prolonged antigen exposure likely contributed to the development of T cell dysfunction 53,54 ; months of induction chemotherapy may have increased tumor cell immunogenicity, eliciting chronic T cell activation with subsequent exhaustion. 64However, we cannot rule out that chemotherapy may have also directly contributed to induction of dysfunctional features in T cells post-treatment. 65Yet, the overall low TOX/TOX2 expression suggests potential retained responsiveness to ICB.
We revealed an abundance of immunomodulatory interactions between T/NK cells and other cells in the TME, which included suppressive macrophages and effector-differentiated Tregs.Effector-differentiation of Tregs, observed in various tumors, has been related to enhanced suppressive capacity, which nominates Tregs as targets for immunotherapeutic interventions. 35,36he identified M2-like-differentiated macrophages have recently moved into the spotlight of immunotherapy: combined anti-GD2/ CD47 therapy to activate M4 resulted in highly promising anti-tumor activity, accompanied by recruitment of M1-like M4 and reduced M2-like M4. 28,29,61 Taken together, the abundance of immunosuppressive cells warrants exploration of immunotherapy combination strategies.
We identified the immunoregulatory NECTIN2-TIGIT axis as a target with high potential for therapeutic intervention.NECTIN2-TIGIT interaction has been reported in other solid tumors, 66,67 underlining its universality in tumor microenvironments.9][70][71] Moreover, combined TIGIT/PD-L1 blockade enhanced antitumor responses in various solid cancers in vitro and in vivo and the observed in vivo responses shared high similarity with ours, resulting in CR in a part of the animals. 68,69,71This heterogeneous response to ICB in syngeneic animal models has been previously related to M4-driven ICB resistance, 72 which again highlights the potential of therapeutic interventions directed at M4 to increase responses to ICB.Moreover, it underscores the importance of identifying biomarkers for treatment response to select eligible patients.Since the patient numbers in our study were too small to perform patient subgroup analyses, this would be a valuable additional analysis in future studies.][75] Our ultimate effort of translation toward an early phase clinical trial, mimicking relapsed/refractory patients with chemotherapyresistant tumors in vivo, demonstrated a survival benefit of adding TIGIT blockade to the currently favored relapse backbone treatment consisting of TEM/IRI+anti-GD2.This paves the way toward clinical development of TIGIT blockade as an alternative ICB in pediatric solid cancers.Moreover, the vast immunoregulatory network in neuroblastoma which we identified provides a rationale and direction for immunotherapy combination strategies.The number and complexity of the immunosuppressive signaling axes offers a likely explanation for the lacking efficacy of mono-immunotherapies. 19,20 To overcome immunosuppressive signaling in T and NK cells, a combination of immunotherapies targeting multiple immunosuppressive pathways simultaneously may be required.For example, the synergistic effect of TIGIT/PD-L1 blockade may be effectuated by their converging restoration of CD226 signaling. 76The already described anti-GD2/ CD47 combination, and a case report of two refractory (I) Detailed immune cell, macrophage, and T cell composition in small and large tumors.*p < 0.05 between small and large, # 0.05<p < 0.1 between small and large.Two-way ANOVA.(J) Immune cell populations as percentage of total live cells (vehicle condition).Mann Whitney U test.(K) TIGIT and PDCD1 expression in single-cell RNA sequencing data (extracted from R2 (r2.amc.nl;Metelitsa-124509-Seurat_cp10k-GSE223071), of isolated CAR-NKT cells (infusion products) pre-and post-infusion into neuroblastoma patients, in responding and non-responding patients 58 .Two-way ANOVA.***p < 0.001, *p < 0.05.Boxplot represents mean, min, and max.Also see Figure S9.neuroblastoma patients treated with anti-GD2/PD-1 combination therapy illustrate the potential of such combinations. 61,77Importantly, future clinical studies investigating combination strategies will have to balance immunotherapy efficacy with the risk of immunotoxicity, and consider biomarkers to predict treatment response.
In conclusion, we have constructed a comprehensive atlas of the neuroblastoma immune environment and identified functionally relevant targets, including the NECTIN2-TIGIT axis, for (combination) immunotherapies.

Lead contact
Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Judith Wienke ( j.wienke-4@prinsesmaximacentrum.nl).

Materials availability
This study did not generate new unique reagents.

Data and code availability
Single-cell RNA-seq data have been deposited at GEO and EGA and are publicly available as of the date of publication.EGA only includes the patient samples, GEO also includes the healthy control data.S1).Information on socioeconomic status was unavailable as this is not regularly recorded in patient files.We did not perform separate analyses based on sex or gender as the relatively small sample size did not allow for this stratification.This might limit the generalizability of our study.The presence of MYCN amplification was determined by fluorescence in situ hybridization at the pathology department.5 adult healthy volunteers donated blood for analysis of healthy immune cells and 8 buffy coats from the blood bank Sanquin were used for in vitro assays.
Studies underlying this paper have received appropriate approval by ethics review boards as per national legislation.Dutch tumor samples were obtained through an institutionally approved research study by the institutional review board of the Erasmus Medical Center Rotterdam (MEC-2016-793), the institutional review board for biobank at the University Medical Center Utrecht (TC18-774) and the biobank committee of the Princess Ma ´xima Center (PMCLAB2020.093,PMCLAB2020.124 and PMCLAB2019.067).The studies were conducted in patients in accordance with the Declaration of Helsinki.Age-appropriate written informed consent from patients and/or parents was obtained prior to inclusion of each patient in the study.

Animals
All experiments were approved by the Institute of Cancer Research Animal Welfare and Ethical Review Body and performed in accordance with the UK Home Office Animals (Scientific Procedures) Act 1986, the United Kingdom National Cancer Research Institute guidelines for the welfare of animals in cancer research and the ARRIVE guidelines. 79 The chemotherapy-resistant model was generated by allograft of neurospheres derived from Th-ALK F1174L /MYCN 129/SvJ transgenic spontaneous tumors into syngeneic animals, which were then subjected to repeat and escalating chemotherapy (VAC -Vincristine, Adriamycin/Doxorubicin, Cyclophosphamide) (detailed in Figures S9A-S9F).The Th-ALKF1174L mutation was introduced as described previously 57 and Th-ALKF1174L founders were derived from CBA3 C57BL/6J mice and genetically crossed with Th-MYCN mice of the 129/SvJ strain.Mice were bred at ICR.For treatment with ICB, the 129 SVJ mouse colony was maintained at UCL ICH Animal Facility.Animals were randomized for treatment using a Latin Square based on tumor volume, as measured by callipers.Mice were all female, and tumors were engrafted between 6-12 weeks of age.
Animals were housed in autoclaved, individually ventilated cages on racks, with a maximum stocking density of 5 animals per cage, 4 for males.Rooms were light cycled, temperature and humidity controlled, and quarterly health screening was performed to maintain a specific-pathogen-free status.Mice were allowed access to sterile food and water ad libitum.For the duration of the treatment with TEM/IRI, mice were given a damp diet in addition to the standard chow, to minimize the likelihood of chemotherapy-driven gastrointestinal toxicity.The damp diet was started on the day before TEM/IRI treatment start and terminated on the final day of TEM/IRI treatment.Animals were then returned to a diet of standard chow only.

Tumoroids
Patient-derived tumoroids (MHC-1 + AMC691T and MHC-1 -AMC691B) were generated as described previously 81 from a female patient.In short, tumor pieces were minced, enzymatically digested with collagenase I, II and IV to achieve a tumor digest.This digest was cultured at 37 C with 5% CO 2 in neuroblastoma tumoroid medium, consisting of Dulbecco's modified Eagle's medium (DMEM) with low glucose and GlutamaxÔ, supplemented with 20% Ham's F-12 Nutrient Mix, B-27 supplement minus vitamin A (503), N-2 supplement (1003), 100 U/mL penicillin, 100 mg/mL streptomycin, 20 ng/mL animal-free recombinant human EGF, 40 ng/mL recombinant human FGF-basic, 200 ng/mL recombinant human IGF-I, 10 ng/mL recombinant human PDGF-AA and 10 ng/mL recombinant human PDGF-BB, until successful sphere formation and expansion of tumor cells.After generation, the tumoroids were frozen in medium with 10% DMSO until further use.For experiments, tumoroids were thawed and cultured in neuroblastoma tumoroid medium at 37 C with 5% CO 2 which was refreshed twice per week.Depending on the growth rate, cells were subcultured 1:2 once or twice per week.Cell authentication was performed regularly by STR profiling.The patient-derived tumoroids AMC691T and AMC691B with GFP-luciferase constructs were generated from these lines previously and viably frozen until further use. 82

Isolation of cells from tumor biopsies
Tumor material was collected in the operating room by Tru-cut biopsy (at diagnosis) or surgical resection (after induction chemotherapy).Preparation of the tumor pieces was started within 4 hours after surgery.The material was minced into pieces <1mm 3 and dissolved by enzymatic digestion with collagenase I, II and IV (2.5 mg/mL) at 37 C with agitation for max. 1 hour, filtered through a 70 mm cell strainer and washed in DMEM.Cells were further separated into a single cell suspension with the NeuroCult dissociation kit according to the manufacturer's protocol (StemcellÔ Technologies, cat#05707), washed and stained with DAPI and DRAQ5.Some samples were additionally stained with antibodies to enrich for T cells (NB124 and NB125 with CD3-PE; NB124 with TCRgd-BV421) or for tumor cells (000GGU, 000GXF, NB106, NB107, NB098, NB125, NB125, NB130, NB152 with anti-GD2).Single live cells were FACS sorted into 384-well plates containing 10 mL of mineral oil, 50 nL of Reverse Transcription primers, deoxynucleotide triphosphates (dNTPs) and synthetic mRNA Spike-Ins on a FACSJazz, FACSAria II or Sony SH800S machine and subsequently spun down, snap-frozen on dry ice, and stored at -80 C until further use.

Isolation of cells from peripheral blood
For single-cell RNA sequencing, peripheral blood mononuclear cells (PBMC) from 5 healthy young adult donors (mean age 28, 2 male, 3 female) were freshly isolated by Ficoll-Paque density centrifugation and stained with fixable viability dye efluor506, CD3-PE and CD56-APC antibodies.Live CD3 + T cells and CD56 + NK cells were FACS sorted into 384-well plates containing 10 ml of mineral oil, 50 nl of RT primers, deoxynucleotide triphosphates (dNTPs) and synthetic mRNA Spike-Ins on a Sony SH800S machine.The plates were spun down, snap-frozen on dry ice, and stored at -80 C until further use.
For killing assays, PBMC from 8 healthy donors (3 for full PBMC, 3 for immune subsets, and 2 for NK cells) were isolated from buffy coats by Ficoll-Paque density centrifugation.The latter cells were either stained with TCRgd-BV421, CD4-FITC, CD8-PE, CD56-APC, and CD3-APC-Cy7 antibodies (for immune subsets) or with CD3-PE and CD56-APC (for NK cells) and FACS sorted on a Sony SH800S machine.

Immunohistochemistry
Immunohistochemistry of human neuroblastoma paraffin-embedded tissue slides with hematoxylin & eosine, anti-CD3 and anti-SOX10 was performed by the pathology department of the PMC.For CD3 and SOX10 staining, tissue slides were pre-treated with a citrate solution for 20 minutes at 100 C, and incubated with the antibodies for 15 minutes at room temperature.The stainings were performed on a BOND immunostainer and visualized with the BOND polymer refine detection kit with a DAB enhancer.The stained slides were subsequently annotated by our neuroblastoma-specialized pathologist.Immunohistochemical staining of mouse neuroblastoma tumors was performed using anti-TIGIT, anti-PD1 and anti-PD-L1 antibodies at WuXi Biologics.
CEL-Seq2 library preparation, sequencing & mapping All samples were processed for total transcriptome amplification, library preparation and sequencing into Illumina sequencing libraries as previously described. 83Paired-end 2x75 bp sequencing read length was used to sequence the prepared libraries using the Illumina NextSeq sequencer.Sharq preprocessing and QC pipeline were applied to process the single-cell RNA-seq data as described. 84Read mapping was done using STAR version 2.6.1 (RRID: SCR_004463) on the hg38 Patch 10 human genome.Function featureCounts (RRID: SCR_012919) of the subread package (version 1.5.2) was used to assign reads based on GENCODE version 26 (RRID: SCR_014966).

CEL-Seq2 quality control
Failed reactions were identified by low levels of ERCC external RNA controls and excluded. 84Furthermore, a liveness threshold was calculated for each plate based on the wells with no cell added, in order to distinguish live cells from dead and/or apoptotic cells. 84his threshold was set to a minimum of 500 transcripts.Genomic:protein-coding read ratio (GPratio) was calculated based on raw counts and cells with a GPratio <20 were removed. 84ext, the percentage of transcripts mapping to the mitochondrial genome was calculated and cells with more mitochondrial-encoded transcripts over nuclear ones were removed.Mitochondrial and ERCC transcripts were removed from the dataset, as well as cells with <1000 nuclear-encoded transcripts or <500 genes.In addition, cells with >150.000nuclear-encoded transcripts were removed.Pseudogenes were removed, as well as genes with low expression -that is either having less than 5 cells expressing the gene or less than two cells with less than two transcripts.A distinct cluster of erythroid lineages was identified based on high levels of hemoglobin complex genes, and was removed.Between-sample variation was minimized by using a standard operating procedure for sample preparation (described above).To further improve cross-sample comparisons, ambient mRNA contamination in individual cells was estimated and removed using DecontX. 85DecontX was run for all samples (batches) individually.Removal of cells with less than 1000 nuclear-encoded transcripts was repeated on the decontaminated counts matrix outputted by DecontX.Afterwards, one cluster still showed a higher contamination score (>0.2), as compared to the other clusters, and was additionally removed.
Single-cell RNA sequencing cluster identification All subsequent analyses were performed using R (version 4.0.2) and the package Seurat (version 3.2.2) with default parameters unless stated otherwise.The SCTransform() function was used to normalize and scale the data, and to identify variable genes.To avoid clustering of cells based on specific cell processes, genes associated with sex (XIST, TSIX, and Y chromosome-specific genes), cell cycle phase, dissociation stress (heat shock proteins; GO:0006986), and activity (ribosomal protein genes; GO:0022626), were removed as described before. 86rincipal component analyses were performed using the filtered lists of variable genes.To study the main cell types in the tumor, the first 35 principal components (PCs) were used to calculate dimensionality reduction using UMAP, and a resolution of 0.5 was used to define clusters using the Louvain method.For immune cell-focused analysis immune cell clusters were subset based on PTPRC gene expression.For in-depth analyses, the respective clusters were subset and UMAP was rerun with an optimal number of PCs and resolution to define subclusters (T/NK: 45 PCs, resolution 0.6; myeloid: 35 PCs and resolution 0.8; mesenchyme: 40 PCs, resolution 0.2).
For the peripheral blood healthy control dataset the first 8 PCs were used for dimensionality reduction, and a resolution of 0.5 was used to define the clusters.For combined analysis of tumor and blood, the two Seurat objects were merged, after which 50 PCs and a resolution of 0.4 were used for dimensionality reduction and defining clusters.For in-depth analysis of the T and NK cells the respective clusters were subset from this combined dataset, UMAP was rerun using 35 PCs and a resolution of 0.5.

Cluster annotation
Cluster annotation was performed with R package SingleR (version 1.2.4), 87 using the HumanPrimaryCellAtlas reference dataset for main cell types, and the NovershternHematopoieticData 88 and MonacoImmuneData 89 reference datasets for immune cell (sub)clusters.Cell annotations were further refined by consulting cluster-specific (up-regulated) differentially expressed marker genes using Seurat's FindAllMarkers function.1][92][93] Malignant and non-malignant cells were distinguished according to three criteria: (1) their inferred CNV profiles (see below); 94 (2) under-expression or absence of different non-malignant cell type marker genes; and (3) high expression of published neuroblastoma-associated genes (Figures 1D and S1B).For tumor cell identification by copy number variation inference the R package inferCNV was run (with default settings except cutoff=0.1,denoise=T and HMM=F) using the immune and endothelial cell clusters as a reference.
To mirror our identified subsets against the subsets of Costa et al. and Verhoeven et al., we used gene signatures from their datasets (see Table S2) and performed gene set enrichment analysis with the differentially expressed genes between subsets in our dataset (Findallmarkers with adjusted parameters, see code).

Differential gene expression analysis
Cluster-specific genes were identified using the Findallmarkers function in Seurat (RRID:SCR_007322) and genes with padj<0.05were considered differentially expressed.Differentially expressed genes between two groups were identified using the Findmarkers function and genes with padj<0.05were considered differentially expressed, as indicated in the figure legends.Volcanoplots were created with the R package EnhancedVolcano.

Cellular composition analysis
For analyses determining the composition of cell types of the individual tumor samples, we only included samples which were sorted in an unbiased manner (DAPI & DRAQ5 for total cells, overall immune cells, and myeloid cells; DAPI & DRAQ5 or anti-CD3 for T cell subsets).For immune cell composition analyses, samples with <10 immune cells were excluded.Barplots were generated using ggplot2 (version 3.3.2).Differences in composition between treatment-naive and treatment-exposed tumor samples were statistically tested using the Wilcoxon rank-sum test.
Pathway and gene set enrichment analysis Pathway enrichment analysis was conducted with the online Reactome portal (https://reactome.org) using differentially expressed genes with padj<0.05.Pathways with >10 identified genes and Bonferroni-corrected P<0.05 were considered statistically significant.
For GSEA, differential expression analysis was performed with Findmarkers for two groups, or Findallmarkers for >2 groups, using adjusted parameters (see code).Genes were pre-ranked by their Fold Change and GSEA was performed using Broad Institute software, by 1000 random permutations of the phenotypic subgroups to establish a null distribution of enrichment score, against which a normalized enrichment scores and multiple testing FDR-corrected q values were calculated. 95Gene sets with an FDR<0.25 were considered significantly enriched, as recommend by Broad Institute.Gene sets were either obtained from MSigDB, from provided data in publications or by analyzing raw data using GEO2R (NCBI tool). 96An overview of used signatures is provided in Table S2.

Interaction analysis
The CellChat algorithm was applied to perform an unbiased ligand-receptor interaction analysis, using the curated ligand-receptor database of CellPhoneDB. 52To identify functionally relevant interactions, we overlaid 1) genes B expressed by cell subset X which were predicted to be membrane-expressed or secreted in the Human Protein Atlas database 97 and had a significant correlation (p<0.05) with the cytotoxicity score or exhaustion score of T/NK cell subset Y and 2) were predicted to have a significant interaction (p<0.05) between cell subset X and T/NK cell subset Y.

Flow cytometry of neuroblastoma tumors
Tumor samples were minced, enzymatically digested as described above, but only with collagenase IV to preserve surface marker expression, and stored as cellular digest in FCS with 10% DMSO in liquid nitrogen until further use.Samples were thawed, washed twice with PBS, stained for 20 minutes at -4 C with Fixable viability dye eFluor 506 (eBioscience) subsequently washed in PBS, stained with the surface antibodies (see key resources table ) for 20 minutes at 4 C, washed in staining buffer (PBS with 2% FCS and 2mM EDTA) and fixed and permeabilized with the Foxp3/Transcription Factor Staining Buffer Set (eBioscience) according to manufacturer's instructions.Cells were then stained with the intracellular and intranuclear antibodies for 30 minutes at 4 C, washed and measured on an Aurora (Cytek) instrument.

In vitro tumor killing assays
Cultured tumoroids were dissociated into single cells with Accutase (Sigma) and mechanical dissociation by pipetting.5000 cells per well were seeded in 100 ml tumoroid medium in a white flat-bottom TC-treated 96-well plate (Corning, cat#3917) and rested for 2-3 days to reform spheres at 37 C and 5% CO 2 .Effector cells from healthy donors, i.e. isolated allogeneic PBMCs, immune subsets, NK cells or PRAME-TCR transduced T cells (the latter were kindly provided by the Nierkens group 98 ), were added to the tumoroids at the indicated effector:target (E:T) ratios in RPMI 1640 (Gibco) supplemented with 10% FBS, 100 U/mL Penicillin/Streptomycin and 2mM% L-Glutamine, with or without anti-TIGIT (clone 10A7) and anti-PD-L1 (clone 6E11; both kindly provided by Roche; 10 ng/mL), recombinant human TGF-b (10 ng/mL), or anti-TGF-b antibody (1 mg/mL).NK cells were primed before the co-culture with IL-2 (1000 IU/mL) and IL-15 (50 ng/mL) overnight.After 6 days (PBMC, immune subsets and PRAME TCR-T cells) or 24 hours (NK cells), supernatants were collected for multiplex immunoassay (for measurement of granzyme B expression by NK cells) 99 and D-luciferin (150 mg/ mL) was added to the culture and incubated for 5 minutes at 37 C.The produced luminescence signal was detected with the FLUOstar Omega microplate reader and normalized against an untreated tumoroid only control to calculate the percentage of tumor killing.

TIGIT and PD-L1 blockade in vivo
In vivo experiments were carried out in Neuro-2a, N1E-115 and N18 syngeneic tumor models.Neuroblastoma cells were routinely subcultured twice weekly.The cells growing in an exponential growth were subcutaneously inoculated at 1x10 6 cells/ml in 0.1 ml PBS for tumor development, and tumors were grown for 7 (N1E-115), 6 (Neuro2a) and 5 (N18) days before start of treatment, until average tumor volume was $100 mm 3 .On day 0 of treatment, body weight and tumor volume were measured, and mice were randomized into study groups with nine animals per arm.Mice were dosed IP with anti-mu-TIGIT (10 mg/kg, first dosing done i.v.) and anti-PDL1 clone 6E11 (10mg/kg), supplied by Roche.Both molecules were administrated 3 times per week for 3 weeks.For TIL analysis tumors were harvested, weighed, and processed to obtain single-cell suspensions; tumor tissues were dissociated using MACS Dissociator (Miltenyi-130-093-235).All samples were analyzed by flow cytometry (BD FACS LSR Fortessa) and >10,000 CD45 + cells were recorded for further analysis.
The chemotherapy-resistant model was generated by allograft of neurospheres derived from Th-ALK F1174L /MYCN 129/SvJ transgenic spontaneous tumors into syngeneic animals, which were then subjected to repeat and escalating chemotherapy (VAC -Vincristine, Adriamycin/Doxorubicin, Cyclophosphamide).These chemotherapy-resistant tumors were subsequently allografted into untreated syngeneic animals which underwent treatment with TEM/IRI (Temozolomide: 3.2 mg/kg IP; Irinotecan: 1.6 mg/kg IP; daily on days 1 through 5:) + anti-GD2 (14G2a, 0.8 mg/kg IP on day 1 and day 5) with or without anti-TIGIT and/or anti-PD-L1 (10 mg/kg IP; twice weekly for three weeks).Mice with tumor volumes <50 mm3 (small tumors) or tumor volumes 50-200 mm3 (large tumors) at the start of treatment (day 10 after engraftment for small tumors and day 14-20 after engraftment for large tumors) were included in the analysis.For TME analysis, mice were culled by humane methods, tumors were harvested and dissected, and tumor pieces were cut by scalpel, washed twice with PBS, and put on a rotator in 5 mL accumax for 10 minutes.Cells were then pressed through a cell strainer and washed twice with PBS.The pellet was resuspended in 1 mL K lysis buffer for 5 minutes and washed twice with PBS.The cells were then counted using a CellCountess and stained according to manufacturer's instructions with antibodies, and run on a BD Symphony cytometer and matched software.

QUANTIFICATION AND STATISTICAL ANALYSIS
Apart from the differential gene expression analysis described above, comparisons between two groups (e.g. for signatures) were made by Mann-Whitney U test.Comparisons between 3 or more groups were analyzed by Kruskal-Wallis with Dunn's post hoc test for multiple comparisons or by 2-way ANOVA for multi-level analyses, as indicated in the figure legends.For heatmap analysis, normalized gene expression was extracted from DotPlot analysis in Seurat, and hierarchical clustering was performed with Ward's method and Euclidian distance.For correlation analyses, Pearson correlation coefficients were computed.FACS data were analyzed using FlowJo V10.8.1 software (LLC).Statistical analysis of FACS data was performed with Graphpad Prism (GraphPad Software).For survival analysis, matched analysis using the log-rank (Mantel-Cox) test was performed.Differences in tumor volumes between treatment groups in in vivo studies were calculated by linear mixed-effects models analysis (see code).

Figure 1 .
Figure 1.The single-cell landscape of neuroblastoma (A) Overview of patients and samples.(B) UMAP of main cell types.(C) Heatmap of top 5 differentially expressed genes per main cell type.(D) UMAP showing expression of neuroblastoma-associated genes.(E) Proportion of different cell types per sample.*Only cells which were sorted with an unbiased FACS sorting strategy were included, which led to exclusion of M761AAA_T, M241AAE_T, and M259AAA_T.''T1'' and ''T2'' refer to paired samples before and after treatment, respectively.(F) Average cellular composition of samples before and after induction chemotherapy.Mixed-effects analysis with Sidak's multiple comparisons test.*p < 0.005.(G) Proportion of mesenchymal and tumor cells in five paired pre-and post-treatment samples.Mann-Whitney U-test.Also see Figure S1.

Figure 3 .
Figure 3. Neuroblastomas are characterized by lymphoid populations with differing degrees of dysfunctionality (A) UMAP and dotplot of T and natural killer (NK) cell subclusters with a selection of their marker genes.(B) Expression of LAG3 and PDCD1 in T/NK cell clusters.*padj < 0.0001 in FindAllMarkers analysis among lymphocytes.(C) Transcription factors associated with effector Treg profile.*padj < 0.0001 in FindAllMarkers analysis among lymphocytes.(D) Volcano plot of differentially expressed genes between the two CD4 T cell populations.(E) Expression of a previously published signature for CD4 T cell dysfunction in melanoma.Mann-Whitney U test. 33(F) Flow cytometric identification of two CD4 T cell populations by expression of PD-1 and IL-7R.(G)GSEA comparing neuroblastoma immune cell clusters (in rows) to previously published gene signatures which identify T cell clusters isolated from either blood, healthy tissues (''normal'') or tumors (in columns).42The color scale indicates for each neuroblastoma cluster the degree of similarity (NES) with published signatures.Also see FigureS3.

Figure 5 .
Figure 5. Immunoregulatory interactions in neuroblastoma(A) Interaction network of main cell types in neuroblastoma constructed with CellChat.52(B) Bubbleplot of most frequent predicted interactions between T/NK subsets and all other cells in the TME.Interactions with >5 partners for NK cells, >2 partners for T cells, and >5 partners for T/NK cells (in at least one subset) are shown.The highest probability per interaction per T/NK cluster is indicated.Interactions shown in Figure5Care highlighted.(C) Bubbleplot of selected, predicted immunosuppressive interactions with T/NK subsets.Interactions with each specific myeloid subset were evaluated and subsequently merged, with the highest probability of each interaction pair depicted in the plot.Also see FigureS6.

Figure 6 .
Figure 6.Targetable immunoregulatory interactions in neuroblastoma(A) Graphical representation of analysis strategy for Figure6B: selection of genes (''B'') expressed by population ''X'' which are involved in a significant ligandreceptor interaction between population ''X'' and T cell subset ''Y'', of which the expression by population ''X'' also significantly correlates with the dysfunction score of T cell subset ''Y''.Genes with at least one significant correlation, with either DUSP4 hi CD4 or CD8 T, were included.(B) Heatmap showing Pearson correlation of expression of genes ''B'' by populations ''X'' with dysfunction score of T cell subsets ''Y''.(C) All predicted NECTIN2-TIGIT interactions with DUSP4 hi CD4 and CD8 T cells in neuroblastoma.(D) Flow cytometric validation of nectin-2 protein expression on neuroblastoma tumor samples.Mean ± SD. (E) Flow cytometric validation of TIGIT protein expression on T cell populations infiltrating neuroblastoma, compared to reference T cells from blood (PB).TPT1 hi CD4 were gated as IL-7R hi PD-1 lo CD4 + cells and DUSP4 hi CD4 were gated as IL-7R lo PD-1 hi CD4 + cells, as shown in Figure3F.Two-way ANOVA with Sidak's posthoc test.Mean ± SD. (F) Correlation of NECTIN2 gene expression with dysfunction score in bulk-RNAseq dataset of SEQC cohort consisting of 498 neuroblastomas (r2.amc.nl;Tumor Neuroblastoma-SEQC-498-RPM-seqcnb1; GSE49710).Also see FigureS7.

Figure 8 .
Figure 8. TIGIT blockade improves survival in a chemotherapy-resistant neuroblastoma model (A) Graphic representation of study setup, adding anti-TIGIT and/or PD-L1 to the standard relapse backbone treatment (Temozolomide/Irinotecan (TEM/IRI) + anti-GD2) in a chemotherapy-resistant, immunologically cold model.

Figure S1 .
Figure S1.The single-cell landscape of neuroblastoma, related to figure 1. (A) UMAP of identified cell clusters, main cell types and sample distribution.Dotted circles highlight the tumor cell clusters.Right panel: number of total cells, tumor cells and immune cells per sample.(B) UMAP with expression of typical neuroblastoma genes (extension of Figure 1D).(C) Copy number profile of the samples, generated using inferCNV, to confirm tumor cell identity.The identified cell clusters (Figure S1A) were used as cell groups, in which immune and endothelial cell clusters were used as reference.*Clusters annotated as tumor cells based on gene expression profile.(D) Volcanoplot of differentially expressed genes (padj<0.05) between tumor cells (malignant) and all other cells (nonmalignant).(E) Violin plot showing the modulescore of HLA class I genes HLA-A, HLA-B, HLA-C, HLA-E, HLA-F and B2M in the main cell types.Kruskall-Wallis + Dunn's.ᵒp<0.0001versus all non-malignant clusters except mesenchyme, which was not signficant.*p<0.0001versus all non-malignant clusters.(F) Correlation between MYCN expression and HLA class I gene score for the three tumor clusters.(G) Cell proportion in MYCN-A versus MYCN-NA samples.(H) UMAP of mesenchyme clusters and sample distribution.(I) Proportion of different mesenchymal clusters per sample."T1" and "T2" refer to paired samples before and after treatment, respectively.*Proportion significantly differs pre vs post treatment over all samples (p<0.05;Mann-Whitney U test).¤Proportion significantly differs pre-vs post-treatment in paired samples (p<0.05;Wilcoxon matched-ranks test).(J) Annotation of mesenchyme clusters with signatures from human protein/cell atlas (HPCA), GSEA using signatures of CAF states and fibroblast types1  , signatures of adrenal cell types 2 , and signatures of follicular dendritic cells (FDC)3  .Distribution pre-vs post-treatment (Tx) per cluster and top differentially expressed gene per cluster.NES = normalized enrichment score.(K) Pathway analysis (Reactome; Bonferroni p-value<0.05) of combined clusters C0-8 of genes which were significantly upregulated (p<0.05)pre-treatment vs posttreatment using all samples.

Figure S2 .
Figure S2.The immune environment of neuroblastoma, related to figure 2. (A) UMAP of sample distribution in immune compartment.(B) UMAP of sample distribution in myeloid cell (left) and lymphoid cell (right) compartments.(C) Flow cytometry gating strategy to identify cell populations in the tumor microenvironment.Merged data from 5 individual tumor samples are shown.(D-E) Gene set enrichment analysis comparing the identified non-immune (D) and myeloid (E) populations in our dataset to identified populations in the datasets of Costa et al. and Verhoeven et al.4,5  .Gene signatures were derived from FindAllMarkers analysis in Seurat (genes with padj<0.05).NES=normalized enrichment score.(F) Heatmap of genes representing antigen presenting/co-stimulatory capacity.Genes are included in module score shown in Figure2C.(G) Heatmap with secreted factors which were 1) specifically upregulated in myeloid cells compared to all other cells (padj<0.05) and 2) among the differentially expressed genes between the myeloid clusters (padj<0.05).

Figure S3 .
Figure S3.Neuroblastoma is characterized by immunosuppressive and dysfunctional lymphoid populations, related to figure 3. (A) Gene set enrichment analysis comparing the identified lymphoid populations in our dataset to identified lymphoid populations in the datasets of Costa et al. and Verhoeven et al. 4,5 .Gene signatures were derived from FindAllMarkers analysis in Seurat (genes with padj<0.05).NES=normalized enrichment score.(B) Immunohistochemical staining of hematoxylin & eosin (H&E), CD3 and SOX10 showing T cell infiltration in a representative neuroblastoma lesion.(C) Expression of effector Treg marker genes which are

Figure S4 .
Figure S4.The immune cell composition and function in neuroblastoma before and after induction chemotherapy, related to figure 4. (A) Composition of overall immune cells per patient before and after induction chemotherapy."T1" and "T2" indicate paired samples before and after treatment.(B) Composition of detailed immune cell clusters per patient before and after induction chemotherapy."T1" and "T2" indicate paired samples before and after treatment.(C) Lymphoid:myeloid ratio before and after chemotherapy.Dotted lines indicate paired samples.Mann-Whitney U test.(D) Myeloid cell percentages before and after induction chemotherapy.Dotted lines indicate paired samples.Mann-Whitney U test.(E) Lymphoid cell percentages before and after induction chemotherapy.Dotted lines indicate paired samples.Mann-Whitney U test.(F) UMAP and dotplot of integrated T/NK subsets from both tumor and reference peripheral blood (PB) of adult healthy donors showing a selection of their differentially expressed and marker genes.(G) Pathway analysis of NK cells with differentially expressed genes pre-versus post-treatment (nominal p<0.05;Reactome, Bonferroni-corrected p-values<0.05).(H) PRF1, GZMB and GNLY expression in two paired samples before and after treatment.The other 3 paired samples did not have sufficient NK cells for a reliable analysis.(I) Gene set enrichment analysis of NK cells in pre-treatment tumors compared to NK cells from post-treatment tumors and blood.NK tumor = pre-and post-treatment NK cells combined.NES=normalized enrichment score.(J) DotPlot of TGF-β1 expression and modulescores of curated signatures for TGF-β downstream signaling in the neuroblastoma TME 6 .(K) Pearson correlation between expression of signature for downstream signaling by the TGF-β receptor (TGFBR) complex (Reactome) in each population shown in figure S4J and the percentage of cells expressing a modulescore of TGF-β receptor 1/2 in those respective populations.(L) Dotplot of NK cell activating and inhibitory receptor expression.(M) Modulescore of NK cell activating and inhibitory receptors shown in Figure S4L.Kruskall-Wallis with Dunn's, ****P<0.0001vs NK PB, ns=not significant.(N) The TIGIT/CD226 and CD96/CD226 gene expression ratio in two paired samples before and after treatment.The other 3 paired samples did not have sufficient NK cells for a reliable analysis.(O) Flow cytometric analysis of TIGIT and CD96 expression in neuroblastoma-infiltrating NK cells compared to blood NK cells.Mann-Whitney U test; ns=not significant.

Figure
Figure S5.T cells in neuroblastoma show increased dysfunctionality post-treatment, related to figure 4. (A) Violin plots of genes associated with cytotoxicity in γδ T and CD8 T cells pre-and post-treatment compared to cytotoxic T cells in reference blood (Tcyt PB).ns=not significant.(B-D) Pathway analysis of differentially expressed genes (padj<0.05) in γδ T (B), CD8 T (C) and CD4 T (D) preversus post-treatment (Reactome, Bonferroni-corrected p-values<0.05).(E) UpsetR 'Venn' diagram of top 50 upregulated pathways pre-treatment compared to post-treatment in the T cell subsets.(F) Gene set enrichment analysis (GSEA) with gene signatures for

Figure S6 .
Figure S6.Immunoregulatory interactions in the neuroblastoma tumor-microenvironment, related to figure 5. (A) Interaction network involving interactions between myeloid cells and T/NK cells constructed with CellChat 9 .(B) Bubbleplot of predicted interactions between myeloid cells and indicated T/NK subsets.Interactions with each specific T/NK subset and each specific myeloid subset were evaluated and subsequently merged, with the highest probability of each interaction pair depicted in the plot (left: myeloid merged; right: T/NK merged) (C) Bubbleplot of predicted interactions between tumor cells and indicated T/NK subsets.(D) Bubbleplot of predicted interactions between mesenchymal cells and indicated T/NK subsets.(E) Interaction network of interactions among Tregs and all other T/NK cell subsets.(F) Bubbleplot involving interactions between Tregs and indicated T/NK cell subsets.

Figure S7 .
Figure S7.The NECTIN2-TIGIT axis is associated with T cell dysfunction, related to figure 6. (A) Dotplot with expression of dysfunction/exhaustion markers CTLA4, LAG3, TIGIT, PDCD1 and HAVCR2 in the indicated T cell subsets.(B) Flow cytometric analysis of CTLA-4, LAG-3, TIGIT, PD-1 and TIM-3 expression on neuroblastoma-infiltrating T cell subsets.TPT1 hi CD4 were gated as IL-7R hi PD-1 lo CD4 + cells and DUSP4 hi CD4 were gated as IL-7R lo PD-1 hi CD4 + cells, as shown in figure 3F.Kruskall Wallis + Dunn's.*P<0.05,**P<0.01,***P<0.001,ns=not significant.(C) Modulescore of the dysfunction/exhaustion genes in the indicated T cell subsets.Kruskall Wallis + Dunn's.****P<0.0001.(D) NECTIN2 'ranking' among genes "B" from figure 6B which correlate positively with dysfunction score of T cells, showing the mean sum of correlations of all interactions involving gene "B" in figure 6B, and the total number of interactions involving gene "B" in figure 6B.(E) Correlation between NECTIN2 expression in IL10hi Mφ and total cells in the TME with exhaustion scores of DUSP4 hi CD4 and CD8 T cells.(F) NECTIN2 gene expression in neuroblastoma TME.*P<0.0001 from Findallmarkers analysis.(G) NECTIN2 gene expression in detailed mesenchyme clusters.*P<0.0001from Findallmarkers analysis.(H) TIGIT gene expressing on T cells infiltrating neuroblastoma compared to T cells from healthy donor blood.*P<0.0001from Findallmarkers analysis.(I) TIGIT and NECTIN2 expression in dataset of Verhoeven et al. 5 (J) TIGIT and NECTIN2 expression in dataset of Costa et al.4

Figure S8 .
Figure S8.Combined TIGIT/PD-L1 blockade enhances immune responses against neuroblastoma in vitro and in vivo, related to figure 7. (A) DotPlot with expression of PD-1 ligands CD274 (PD-L1) and PDCD1LG2 (PD-L2) in neuroblastoma.(B) Flow cytometric analysis of PD-L1 expression in the neuroblastoma tumor-microenvironment. (C) Gene set enrichment analysis of selected signatures for PD-1 and/or TIGIT downstream signaling.KO=knock-out, GC-Tfh=germinal center follicular helper T cells.NES=normalized enrichment score.(D) Gating strategy for flow cytometric expression analysis of Nectin-2 and PD-L1 expression on GFP-positive tumoroids as shown in Figure 7C.(E) Flow cytometric analysis of PD-1 and TIGIT expression on healthy donor PBMC cocultured with neuroblastoma tumoroids for 6 days.'PBMC w/o tumor' was measured at timepoint 0 before coculture.(F) Percentage of tumoroid (AMC691T) killing (=100-normalized luminescence) in vitro after 6 days of co-culture with PRAME-TCR transduced T cells at E:T ratio 1:300.Luminescence values were normalized against condition with untreated tumoroids only.n=4-5 technical replicates.Kruskall-Wallis with Dunn's.*p<0.05;ns=not significant.(G) Percentage of tumoroid (AMC691T) killing (=100normalized luminescence) in vitro after 6 days of co-culture with FACS-sorted immune subsets from healthy donor peripheral blood at E:T ratio 1:3.Luminescence values were normalized against condition with untreated tumoroids only.n=3 donors.Two-way ANOVA with Sidak.*p<0.05;ns=not significant.w/o=without.(H) Tumor volumes at treatment start per condition for the three in vivo models.(I) Immunohistochemical staining of TIGIT, PD-1 and PD-L1 in the three in vivo models (vehicle condition).(J) Body weight of N1E-115, Neuro2a and N18 mouse models (n=6 per group) treated with anti-TIGIT and/or anti-PD-L1 up to day 15 of treatment.(K-N) Flow cytometric analysis of the TME in vivo in n=3 mice per treatment condition, treated with anti-TIGIT and/or anti-PD-L1 for 7 days.(K) Percentage of CD45+ immune cells of total live cells.(L) Total immune cell composition across treatments.Models were pooled.*P<0.05versus combination, ****P<0.0001versus combination.2-Way ANOVA.(M) T cell composition across treatments.Models were pooled.*P<0.05versus combination.2-Way ANOVA.(N) Macrophage composition across treatments.Models were pooled.2-Way ANOVA.

Figure S9 .
Figure S9.TIGIT blockade improves survival in a chemotherapy-resistant neuroblastoma model, related to figure 8. (A) Generation of chemotherapy-resistant neuroblastoma model from Th-ALK F1174L /MYCN 129/SvJ models repeatedly treated with VAC (Vincristine, Adriamycin/Doxorubicin, Cyclophosphamide).(B) Summary of dose escalation treatment schedule to achieve resistance.Red arrows indicate when single dose VAC was administered.The number of days taken for the tumor to relapse is indicated below each cycle.At the end of each cycle, the allograft tumor was re-passaged into another animal before treatment started.(C) Percentage relative tumor volume change 7 days after treatment with single dose of VAC at 32 mg/kg or 80 mg/kg.(D) Time taken for allograft tumor to relapse after treatment with single high dose of VAC (80 mg/kg VAC) in Cycle 5 and Cycle 6. Tumor cells after Cycle 6 were to establish allografts for therapeutic studies.(E) Macroscopic confirmation of tumor presence for small tumors.(F) Histological confirmation of tumor cell presence for small tumors.(G) TSNE representation of flow cytometric analysis of PD-L1 and Nectin-2 expression in the TME of small tumors (vehicle condition).(H) Body weight of mice with small tumors (n=7 per group) treated with TEM/IRI, anti-GD2 and/or anti-TIGIT and/or anti-PD-L1 up to day 80 of treatment.(I) Tumor volumes at start of treatment per treatment condition for mice with large and small tumors.(J) Survival analysis comparing response to TEM/IRI (T/I) +anti-GD2 (aGD2) in small and large tumors.(K) TSNE representation of flow cytometric analysis of Nectin-2 expression in the TME of large tumors (vehicle condition).(L) Flow cytometric analysis of immune cell subsets as a fraction of total live cells in the tumor.Mann Whitney U test.N=4 animals per condition.
Accession numbers are listed in the key resources table.All original code has been deposited at Bitbucket and is publicly available as of the date of publication.DOIs are listed in the key resources table.Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.From the Princess Ma ´xima Center in Utrecht, The Netherlands, 19 patients with neuroblastoma were included between September 2017 and November 2020.Samples were obtained from patients during diagnostic biopsy pre-treatment (n=10) or during surgical resection after ca. 5 months of chemotherapeutic treatment (n=14).Demographic (e.g.sex, age, race), clinical (e.g.treatment) and histological information, including tumor staging by INRG, INRGSS and INSS was collected (Table ,80For in vivo experiments with syngeneic mouse tumor models Neuro-2a, N1E-115 and N18, approximately six to eight-week-old female A/Jax-000646 mice were ordered from The Jackson Laboratory.The animals had a starting body weight of 21±1.7 (mean±SD) gram.