Clinically Relevant Cytotoxic Immune Cell Signatures and Clonal Expansion of T-Cell Receptors in High-Risk MYCN-Not-Amplified Human Neuroblastoma

Neuroblastomas are pediatric solid tumors originating from neurons of the sympathetic nervous system. The disease course can vary drastically, ranging from relentless progression of high-risk tumors despite aggressive treatment to spontaneous regression of low-risk stage 4S tumors in infants with no cytotoxic therapy interventions (1, 2). Despite recent FDA approval of anti–GD2 antibody-based immunotherapy for high-risk neuroblastoma patients in first remission, the cure rates remain less than 50% accompanied by significant morbidity in survivors (3). Previous large-scale sequencing by us and others has revealed a low somatic mutation burden in this disease with only a few recurrently mutated genes (4–6), imposing a substantial challenge for developing new targeted therapies for a large proportion of patients with high-risk neuroblastoma. Therefore, there remains a critical need to identify effective novel therapies for these aggressive enigmatic tumors.

Immune checkpoint inhibitor therapies have recently shown significant antitumor activity in several adult human cancers (7–9). Their activity is thought to be correlated with higher tumor mutation burden presumably because of the greater probability of infiltrating tumor-reactive T cells recognizing neoantigens presented in these tumors (8). Because neuroblastomas typically have a low somatic mutation burden at the time of diagnosis, it is postulated that there will be fewer T cells in the tumor microenvironment (6). However, recent studies by us and other groups reported that intensive chemotherapy and radiotherapy dramatically increase mutation burden in the relapsed setting (10, 11). Furthermore, DNA-damaging agents in conjunction of dinutuximab have been reported to show a survival advantage for patients with relapsed/refractory neuroblastoma, suggesting that immunotherapy may be further exploited for future therapies (12). Remarkably, in synovial sarcoma, a fusion gene–driven cancer with a low mutation burden, there have been reports of significant responses to immunotherapy using patients' own T cells engineered to express a T-cell receptor (TCR) against NY-ESO-1, a cancer testis antigen expressed on many of these sarcomas (13). These observations and several other studies indicate that immunotherapies directed against tumor-associated antigens other than mutation-induced neoantigens may be a powerful therapeutic modality, even in cancers with low mutation burden (13–15). We therefore performed deep RNA sequencing (RNA-seq) in a cohort of neuroblastomas to characterize the tumor-infiltrating immune cells and determine if their signatures are associated with clinical subgroups, molecular characteristics, and survival outcomes. An independent previously published neuroblastoma RNA-seq dataset (16) was used to validate our findings.

Informed consent from each participant or guardian was obtained by qualified investigators at Children's Oncology Group institutions. Patient's identification was anonymized, and the study was conducted according to the U.S. Common Rule and after approval by the local Institution Review Boards. As an integral part of the Therapeutically Applicable Research to Generate Effective Treatment (TARGET) initiative to genomically profile high-risk neuroblastoma (5, 17), we sequenced the mRNA of 150 primary neuroblastoma tumors obtained prior to any therapy (Table 1 and Supplementary Table S1). Of these, 139 tumors have matched DNA sequencing data from whole-exome and/or whole-genome sequencing (Supplementary Fig. S1). In addition, we used previously published RNA-seq data for 190 patients with neuroblastoma (obtained prior to therapy) as an independent cohort to validate our findings (ref.16; validation cohort; Table 1 and Supplementary Table S1).

Only RNAs with RNA Integrity Number greater than 6.0 were used in this study. RNA-seq libraries were prepared using Illumina TruSeqRNA Sample Preparation V2 kits according to the manufacturer's protocol (Illumina), and 100 bp paired-end sequencing was performed using Illumina HiSeq 2000 sequencers. Each sample was processed through a RNA-seq data analysis pipeline where reads were mapped to the ENSEMBL human genome GRCh37 build 71 using TopHat2 (18). Cufflinks was then used to quantify gene expression as Fragments Per Kilobase of transcript per Million mapped reads (FPKM; ref. 19). The final measure of gene expression level is log₂(FPKM+1). FASTQ files for RNA-seq data are available at the NCI's Genomic Data Commons (https://portal.gdc.cancer.gov) and the high-level data at https://ocg.cancer.gov/programs/target/data-matrix.

The independent validation RNA-seq data used in this study were derived from a published RNA-seq dataset generated using the Illumina HiSeq 2000 platform (16). Gene expression levels in log₂(FPKM+1) and the associated clinical information for these 190 patients (Table 1) were downloaded from the Gene Expression Omnibus database entry GSE49711.

We used consensus clustering for genes with nonzero expression variability in the dataset using the R package ConsensusClusterPlus for class discovery based on the comparison of gene expression profiles (20). Genes located on the Y-chromosome were also removed in order to minimize gender-related bias. Expression level of each gene was then median-centered, and the similarity between expression profiles was quantified using Pearson correlation. Clustering was performed using the average hierarchical clustering (the package's default options), and consensus clustering was run for 10,000 iterations with 80% of the genes and 80% of the samples randomly selected on each iteration. The empirical cumulative distribution function was constructed for a range from 2 to 10 consensus clusters (Supplementary Fig. S2).

The expressed somatic mutations were called using the following approach. First, we used GATK to call RNA variants using the best practices workflow (https://software.broadinstitute.org/gatk/guide/article?id=3891). To ensure the quality of RNA variants, we only kept variants that had a quality score of ≥ 50 with total number of reads ≥ 10, the number of variant reads ≥ 3, and variant allelic fraction ≥ 0.15. To further reduce potential false-positive calls, the expressed somatic mutations were then filtered using nonsynonymous variants that were called in its corresponding tumor DNA sequenced by whole-genome (https://ocg.cancer.gov/programs/target/data-matrix) or whole-exome sequencing (5). Therefore, samples that have only RNA-seq but no DNA sequencing data were excluded from this analysis. A somatic mutation is considered to be called by whole-genome sequencing if it was assigned a numeric somatic score. For whole-exome sequencing, a somatic mutation was required to have at least 3 variant reads and ≥ 0.15 variant allele fraction in tumor DNA, ≤ 1 variant read in normal DNA, and ≥10 total reads in both normal and tumor DNA. One hot-spot mutation in the MYCN gene (MYCN:NM_001293228:exon2:c.C131T:p.P44L) from sample PASUML was added from RNA-seq data manually. The final list of expressed somatic mutations contains 703 mutations in 132 samples (Supplementary Table S2).

Single-sample gene set enrichment analysis (ssGSEA) was performed to calculate gene set enrichment scores for each sample (21) using the ssGSEAProjection module from the GenePattern server (https://genepattern.broadinstitute.org/gp; ref. 22). The MYCN-signature gene set consists of 157 genes up- or downregulated by MYCN (23). The immune score was calculated using the ESTIMATE method (24), and 22 immune cell–specific gene sets were derived from a published study (25). We also used a cytolytic signature, an average expression of six genes, 5 granzymes (GZMA, GZMB, GZMH, GZMK, and GZMM) and perforin (PRF1), similar to what was described by Rooney and colleagues (26), to quantify the level of cytotoxic immune cell activity. Genes that did not have approved HUGO Gene Nomenclature Committee symbols were excluded from the ssGSEA analysis. The ssGSEA scores for TARGET and the validation samples were calculated using a subset of genes shared by both datasets to facilitate the direct comparison of the scores across the datasets (Supplementary Table S3).

TCR sequencing was performed using SMARTer Human TCR a/b Profiling Kits according to the manufacturer's instruction (Takara Bio USA). First, full-length cDNA was generated by reverse transcription from total RNA using an oligo dT primer and a SMART-Seq v4 oligo (template switching oligo). Then, TCR targets were amplified by PCR with a SMART primer and a TCR constant region primer (a and/or b). Final TCR sequencing libraries were generated by a second round PCR using Illumina sequencing oligos and indices. TCR libraries were sequenced on Illumina MiSeq sequencers to generate 2 × 300 base pairs paired-end reads (Illumina). FASTQ files from TCR sequencing or RNA-seq experiments were used to extract TCR sequences and calculate TCR counts using MiXCR 2.1.6 as described (27). TCR beta-chain sequences were used in the analyses to identify the TCR clones.

Previously constructed tissue microarrays (TMA) of neuroblastomas (28) were evaluated for CD8 expression by immunohistochemistry. The TMA was constructed with duplicate 1.00 mm cores from each patient. Immunohistochemistry for CD8 was performed with an antibody purchased from Leica Biosystems (Clone 4B11). CD8 expression was scored as (0) none, (1) rare scatted CD8⁺ cells (less than 5), (2) occasional CD8⁺ cells (less than 10), (3) aggregates of CD8⁺ cells, 10 or more cells with only a single aggregate of 2 or more CD8⁺ cells within 10 μm proximity to each other, and (4) multiple aggregates of CD8⁺ cells. The average of the two independent cores was taken and rounded up. Two-sided Wilcoxon Rank-Sum test was used to compare MYCN-not-amplified (MYCN-NA) and MYCN-amplified (MYCN-A) tumors.

Detailed survival analysis can be found in the Supplementary Materials.

RNA-seq was performed on 150 neuroblastoma tumors as part of the TARGET initiative. A median of 158 million reads per sample was obtained, and 90% of samples yielded >100 million reads (Supplementary Fig. S2). To identify clinically meaningful signatures in this cohort, we first took an unbiased approach using consensus clustering of the gene expression profiles. Stage 4S neuroblastomas (n = 21), which occur exclusively in patients <18 months of age at diagnosis and have reliably undergone complete spontaneous regression of disseminated disease without intervention, were used as a comparator group for the high-risk tumors (n = 129). The optimal solution was 4 clusters, because the shape of the CDF curves did not change much beyond this number (Supplementary Fig. S3). The cluster membership of each of the four groups was associated with distinct molecular and clinical characteristics and survival probabilities (Fig. 1A and B, P < 0.01). MYCN gene amplification status, patient's age, and stage 4S had the largest impact on the composition of the clusters. For instance, cluster 2 contained mostly stage 4S or younger patients (<18 months) with the best patient outcome, and cluster 1 mainly consisted of MYCN-A tumors with significantly poor outcome, in keeping with clinical observations. Clusters 3 and 4 contained the majority of MYCN-NA tumors, and both had poor overall survival like those with MYCN-A tumors. We next explored if there were any differences in infiltrating immune signatures in these 4 clusters using the immune score described by Yoshihara and colleagues (ref. 24; Supplementary Table S3). Cluster 1 which contained most MYCN-A tumors had significantly lower immune scores compared with clusters 2, 3, or 4 (adjusted P < 0.001; Supplementary Table S4) as recently described in an analysis of these data by others (29). However, we found that both clusters 3 and 4 have significantly higher immune scores than clusters 1 and 2, but no significant differences between clusters 3 and 4 (adjusted P = 0.89; Supplementary Table S4).

About 60% of high-risk neuroblastoma are MYCN-NA tumors. Valentijn and colleagues have reported an increased MYCN/MYC activity in a subset of MYCN-NA neuroblastomas, and the activity is associated with poor patient outcomes (23). Therefore, we next attempted to confirm that the reported MYCN-signature consisting of 157 genes was associated with patient outcomes in our cohort. Consistent with their study, we found that this gene signature was significantly associated with the inferior outcomes for high-risk patients with MYCN-NA tumors both in the TARGET cohort (P < 0.0001; Fig. 1C) and in an independent validation cohort (ref. 16; P < 0.05; Fig. 1D). Age of the patients could not explain this association as there was no significant age difference between the MYCN-signature high and low groups in both cohorts (Supplementary Fig. S4).

We next explored tumor-infiltrating lymphocyte (TIL) signatures in high-risk MYCN-NA tumors to determine if there were differences between the MYCN-signature high and low tumors, comparing with those of MYCN-A or 4S. Interestingly, both MYCN-signature high and low groups of MYCN-NA tumors had significantly higher immune scores than MYCN-A tumors in both TARGET and validation cohorts (adjusted P < 0.01, Fig. 2, top plots), indicating elevated TILs in MYCN-NA tumors compared with MYCN-A tumors.

To investigate which immune subtypes were present in these high-risk neuroblastomas, we performed ssGSEA (21) using the CIBERSORT immune cell–specific gene sets (ref. 25; Supplementary Table S3). In keeping with the immune score, MYCN-NA tumors had significantly higher TIL subset scores, in particular activated natural killer (NK) and CD8⁺ T cells compared with MYCN-A tumors (adjusted P < 0.01) with no significant difference between the MYCN-NA tumors with low and high MYCN-signature in both patient cohorts (Fig. 2). We next examined the RNA-seq data for evidence of immune cell–mediated cytolytic activity by quantifying the average gene expression level of perforin (PRF1) and 5 granzymes (GZMA, GZMB, GZMH, GZMK, and GZMM) using a similar approach as Rooney and colleagues (26). Consistent with the immune, NK, and CD8⁺ T-cell scores, the cytolytic scores for MYCN-NA tumors were significantly higher than those of MYCN-A tumors (adjusted P < 0.02), and there was no significant difference between the MYCN-NA tumors with low and high MYCN-signature in both cohorts (Fig. 2). Lower NK, CD8⁺ T-cell, and cytolytic scores were also observed in 4S tumors. These results indicate elevated functional cytotoxic TILs in the MYCN-NA high-risk neuroblastomas compared with the MYCN-A and 4S tumors.

We next examined if the high immune infiltration in tumors was associated with the number of expressed somatic mutations. Interestingly, compared with MYCN-A tumors, MYCN-NA tumors with a high MYCN-signature had increased expressed nonsynonymous mutations; however, MYCN-NA tumors with a low MYCN-signature had the lowest mutation burdens despite the highest immune scores among these three subtypes (Supplementary Fig. S5 and Supplementary Table S2). Therefore, there was no clear correlation between expressed somatic mutation burden and immune signatures. Unfortunately, the mutation data were not available in the validation cohort for the comparable analyses of possible association of mutation load and immune score.

MHC I complex on the surface of all nucleated cells is required to present endogenous cellular antigens to circulating T cells, and the expression of peptide-MHC I (pMHC I) has been reported to be lost or downregulated in many tumors including neuroblastoma to evade host immune surveillance (30, 31). Consistent with these previous reports, the expression of HLA-A, -B, and -C genes was significantly reduced in MYCN-A neuroblastoma in both the TARGET and validation cohorts (Fig. 2, bottom plots). In addition, we confirmed that induction of MYCN expression in a MYCN-NA neuroblastoma cell line resulted in a significant suppression of MHC l genes (Supplementary Fig. S6; ref. 32). Therefore, these data suggest that MYCN-A neuroblastoma may be less immunogenic to T-cell immune surveillance through downregulation of MHC I complex.

The gene expression profiles of infiltrating immune cells in tumor samples have been shown to be prognostic in many human cancers (33). We therefore developed a method, termed K-M optimization procedure (see Supplementary Methods), to identify optimal thresholds of immune signatures associated with outcome for patients with MYCN-NA tumors. By this method, 6 of the 22 immune signatures were significantly associated with patient outcome for tumors with a high MYCN-signature, but none in those with a low MYCN-signature (Fig. 3A). Notably, the activated NK-cell signature along with CD8⁺ T-cell signature exhibited the most significant association with outcome for patients having tumors expressing a high MYCN-signature (Fig. 3A). Out of the 22 immune cell signatures, the cytolytic signature was highly correlated with NK-, T-, and B-cell signatures in the TARGET samples, but only moderately correlated with those of mast cells or macrophages (Fig. 3B). As with the activated NK-cell signature, the cytolytic signature was also significantly associated with outcome of patients with MYCN-signature high tumors in both TARGET and validation cohorts (Fig. 3C and D). Therefore, these findings strongly suggest that activated infiltrating cytotoxic immune cells affect patient outcome in high-risk patients with MYCN-NA tumors expressing a high MYCN-signature.

Because of the important roles of cytotoxic T cells in tumor surveillance and killing, we next examined if we could detect the expression of TCR from our deep RNA-seq data and validated our findings using directly TCR sequencing from two tumors (PAPUEB and PASWVY; see Materials and Methods). We found that the RNA-seq could detect the same highly expressed TCRs as determined by direct TCR sequencing (Supplementary Fig. S7), demonstrating that the RNA-seq can identify the more abundant T-cell repertoires in tumor total RNAs.

We therefore examined if the TCRs detected by the RNA-seq in the tumors were correlated with the immune and the CD8⁺ T-cell scores. We found that the number of total unique TCR clones was significantly correlated with immune or CD8⁺ T-cell scores across all the tumor samples (Fig. 4A). Consistent with the immune and CD8 T-cell scores, the TCR total clone counts were significantly higher in MYCN-NA neuroblastomas than those in MYCN-A tumors (P < 0.0001, Fig. 4B). In addition, the TCR total clone count was also significantly associated with outcome for high-risk patients with MYCN-NA neuroblastoma expressing a high MYCN-signature (P < 0.05; Fig. 4C). Interestingly, we observed a small group of neuroblastomas having clonal expansion (higher TCR clone percentage per sample, e.g. >1%; and higher TCR clone expression, e.g. total TCR counts per clone >10), and they were all MYCN-NA tumors (blue box in Fig. 4D). Our data thus show that MYCN-NA tumors have high TILs, TCR clone counts, and T-cell clonal expansion, and higher TCR clone counts are associated with a moderate but significant survival advantage. Because the sequence data were not available for the validation cohort, we could not confirm these findings in that dataset.

To further validate the involvement of infiltrating CD8⁺ T cells in MYCN-NA neuroblastoma, we performed immunohistochemistry for CD8 on an independent neuroblastoma TMA described previously (28). A total of 80 neuroblastoma samples of all risk patients were examined using a qualitative metric of number of CD8⁺ cells (Supplementary Table S5; see Methods). MYCN-NA neuroblastomas had significantly higher T-cell infiltration than MYCN-A tumors (P < 0.001, Fig. 4E and F). This was confirmed in high-risk–only samples (n = 39), with higher CD8 scores in MYCN-NA neuroblastomas compared with MYCN-A tumors (P < 0.05, Supplementary Fig. S8). Therefore, the results of immunohistochemistry are confirmatory to our findings with the RNA-seq data, supporting that there are significantly elevated TILs in MYCN-NA neuroblastomas.

Exhaustion of CD8⁺ T and NK cells has been recently identified as an important mechanism by which tumors escape host immunity (34, 35). We examined a panel of 7 exhaustion marker genes (LAG3, CD160, CD244, HAVCR2, CTLA4, TIGIT, and PDCD1) together with genes in the cytotoxic immune cell signatures including activated NK, CD8⁺ T cells, and the cytolytic signature in MYCN-NA tumors expressing a high MYCN-signature. These immune cell genes have increased expression in tumors with high cytotoxic cell signatures and better outcome (Fig. 3C and D, and 5A), signifying the importance of functional infiltrating cytotoxic cells in tumor progression. Specially, we observed that 5 of 7 of these exhaustion marker genes were expressed at significantly higher levels in tumors with higher cytotoxic activities in both the TARGET and validation datasets (Fig. 5B), indicating an elevated level of immune exhaustion in conjunction with increased cytotoxic immune activities in these tumor samples. Consistently, these 5 exhaustion marker genes were also expressed at higher levels in MYCN-NA tumors than MYCN-A tumors which had lower cytotoxic immune signatures (Supplementary Fig. S9). Among them, the therapeutically targetable immune inhibitory receptor genes CTLA4 and PDCD1 (PD1) were significantly upregulated in the tumors with higher cytolytic activities.

TILs are indicative of host immune response to tumors, and they have been reported to be predictive of clinical outcomes for patients with cancers including neuroblastoma (26, 29, 33, 36). We investigated TILs in high-risk neuroblastoma using deep RNA-seq on tumors which surveys gene expression of different cell types including tumor, infiltrating immune, and other stromal cells. Complementing a previous report of the roles of tumor-associated macrophages in metastatic neuroblastoma (36), we have shown an elevated activity of cytotoxic immune cells including NK and CD8⁺ T lymphocytes in MYCN-NA neuroblastomas, which account for approximately 60% of all high-risk patients, compared with MYCN-A and 4S tumors. Mutation burden has been reported to be predictive of response to immune checkpoint blockade in human cancers (8, 37, 38). Interestingly, we did not find a direct association between the number of expressed somatic mutations and the cytotoxic immune signatures. Instead, we observed a significant correlation between the expression of genes encoding for MHC I and immune signatures including that of CD8⁺ T cells. Although we observed a slightly higher expressed somatic mutation burden in the MYCN-NA tumors with a high MYCN-signature, their mutation burden remains substantially lower than adult cancers with high mutation burden that exhibit a greater likelihood to respond to immune checkpoint blockade (8, 37, 38). Thus, it is perhaps not surprising that we did not see a clear association between the mutation burden and tumor-infiltrating immune cell activity in all high-risk neuroblastomas in this study. On the contrary, MHCs display endogenous antigens on almost all human nucleated cell surfaces enabling T cells to recognize foreign antigens (resulting from infection or mutated genes) through interaction between peptide-MHC and TCRs (30). MYCN and its homolog MYC have been reported to downregulate HLA genes in neuroblastoma and melanoma respectively (31, 39), and this finding was also confirmed by our study. The loss or downregulation of HLA gene expression is commonly seen in many human cancers as one mechanism of escape from host immune surveillance (30). It is consistent with the observation that MYCN-A neuroblastomas in general have a lower level of tumor-infiltrating T cells as reported in this study and by others (29, 40). Furthermore, a global suppression of IFN and proinflammatory pathways by MYCN also contributes to a T-cell–poor microenvironment in MYCN-A neuroblastomas (40). These observations suggest that MYCN-A neuroblastomas are intrinsically less immunogenic and are unlikely to respond to immune checkpoint inhibitors without additional therapies to induce the expression of MHC class I proteins.

In addition to the discovery of elevated TILs in the MYCN-NA neuroblastomas, we also found that the presence of cytolytic immune cell signatures in tumor samples obtained at diagnosis was associated with improved overall survival in a subgroup of high-risk patients with MYCN-NA tumors expressing a high MYCN-signature. However, these same immune signatures were not predictive for patients with MYCN-NA tumors expressing a low MYCN-signature. This may be due to the already better survival for this group of patients (23), and they often had higher TILs in their tumors as found in this study. One of the exciting observations of this study is the presence of infiltrating CD8⁺ T cells in MYCN-NA tumors, supported by CD8⁺ T-cell signatures, TCR counts, and immunohistochemistry staining of CD8 on a neuroblastoma tissue array. Intriguingly, analyses of TCR sequences in the RNA-seq data demonstrated a potential clonal expansion of T cells in some neuroblastomas, preferentially in the MYCN-NA tumors. Although at this time we do not have definitive evidence for the specificity of these expanded T-cell clones against neuroblastoma cells, we postulate that the clonal expansion is the consequence of T cells recognizing tumor-associated antigens (such as GD2, mutation-induced neoantigens, cancer/testis antigens, or other neuroblastoma-specific developmental genes; ref. 41). This warrants future TCR sequencing of tumors and functional validation of expanded clones to develop these for adoptive cell therapy. In conjunction with CD8⁺ T-cell infiltration, we found a higher expression of cytolytic effector genes such as perforin (PRF1) and granzymes (GZMs), suggesting that these TILs exhibited cytolytic activity. These results suggest that the presence of functional and active cytotoxic immune cells may contribute to why patients with MYCN-NA tumors expressing a high MYCN-signature and elevated cytotoxic immune cell signatures have a moderately but significantly improved outcome over those with lower immune cell signatures. This also suggests that, even at the time of diagnosis, the presence of activated tumor-infiltrating cytotoxic immune cells may play a role in the clinical biology of this subtype of neuroblastomas.

Exhaustion of cytotoxic immune cells is another important mechanism of tumor escape from host immune surveillance, and it is characterized by the inability of cytotoxic T cells to be activated, induce cytokines, proliferate, and result in cytolysis of target tumor cells (42). Exhaustion is accompanied by the expression of inhibitory cell surface receptors, which are effective targets for immune checkpoint blockade therapies to reactivate these cytotoxic immune cells (42, 43). Anergy is another potential cause of unresponsiveness of T cells against NB tumor cells and is associated with high CTLA4 expression (44). Anergy of T cell in tumors is not as well understood due in part to a less well-defined set of markers selectively expressed in anergic cells (42, 45). However, here we observed elevated expression of several other markers not typically associated with anergy including HAVCR2 (TIM3), CD160, and TIGIT. Importantly, the upregulation of these genes was associated with a higher cytolytic signature in MYCN-NA high-risk patients with a high MYCN-signature that had better outcome. We have recently confirmed by immunohistochemistry that the expression of PD-L1, the ligand for PD1, is associated with poor survival for neuroblastoma patients (46). This may explain why even though this group of patients had a modest prolonged survival, many of them eventually succumbed to disease. In these cases, it is also possible that high MYC activity reflected by the MYCN-signature may directly result in a suppressive tumor microenvironment (40). These observations lead to an intriguing hypothesis that these inhibitory receptors may prevent immune cells from maximal killing of tumor cells, thereby accounting for patient death due to cancer despite the elevated immune signatures. If so, this subset of patients with both cytotoxic immune signatures and augmented expression of inhibitory receptors may derive significant potential benefit from immune checkpoint inhibition.

In conclusion, our findings of the association of improved survival with increased signatures for activated NK cells, CD8⁺ T cells, cytolytic activity, clonal expansion of TCRs, and exhaustion markers suggest that newly diagnosed high-risk neuroblastoma patients with MYCN-NA tumors may benefit from multiple modalities of immunotherapy incorporated into their front-line treatment.

No potential conflicts of interest were disclosed.

The content of this publication does not necessarily reflect the views or policies of the Department of Health and Human Services, nor does mention of trade names, commercial products, or organizations imply endorsement by the U.S. government.

Conception and design: J.S. Wei, I.B. Kuznetsov, S. Zhang, J.M. Maris, J. Khan

Development of methodology: J.S. Wei, I.B. Kuznetsov, S. Zhang, D. Catchpoole, S.M. Hewitt, J. Khan

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J.S. Wei, Y.K. Song, S. Asgharzadeh, A. Yuksel, D. Catchpoole, S.M. Hewitt, R. Seeger, J. Khan

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J.S. Wei, I.B. Kuznetsov, S. Zhang, Y.K. Song, S. Sindiri, X. Wen, R. Patidar, S. Najaraj, A. Walton, A. Yuksel, S.M. Hewitt, J.M. Maris, J. Khan

Writing, review, and/or revision of the manuscript: J.S. Wei, I.B. Kuznetsov, S. Zhang, S. Asgharzadeh, S. Sindiri, J.M. Guidry Auvil, D.S. Gerhard, A. Yuksel, S.M. Hewitt, P.M. Sondel, R. Seeger, J.M. Maris, J. Khan

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J.S. Wei, S. Sindiri, X. Wen, J.M. Guidry Auvil, J. Khan

Study supervision: J.S. Wei, D.S. Gerhard, J.M. Maris, J. Khan

The authors thank the Children's Oncology Group for providing samples and related clinical information for this study and all TARGET investigators for support of this study. They also thank Dr. Julie Park, Hongling Liao, Dr. Berkley Gryder, Dr. Michael Krivanek, Dr. Bob Hawley, and the Center for Cancer Research Sequencing Facility (CCR-SF) for their technical and scientific support. This study utilized the high-performance computational capabilities of the Biowulf Linux cluster at the NIH (http://biowulf.nih.gov), and was supported by the Intramural Research Program of the NIH, National Cancer Institute, Center for Cancer Research, and NIH grant RC1MD004418 to the TARGET consortium, and The Stand Up to Cancer - St. Baldrick's Pediatric Dream Team Translational Research Grant (SU2CAACR-DT1113). Stand Up To Cancer is a division of the Entertainment Industry Foundation. Research grants are administered by the American Association for Cancer Research, the Scientific Partner of SU2C.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.