Verication of Genetic Differences and Immune Cell Inltration Subtypes in Neuroblastoma Tumor Microenvironment for Immunotherapy

Background: The tumor microenvironment (TME) has achieved remarkable results in the research of cancer progression in the past few years. it is crucial to understand the nature and function of TME in tumors because of precise treatment strategies for individual cancers having received widespread attention, including immunotherapy. The immune inltrative proles of neuroblastoma (NB) have not yet been completely illustrated. The purpose of this research is to analyses tumor immune cell inltration (ICI) in the microenvironment of NB. Methods: We applied CIBERSORT and ESTIMATE algorithms to evaluate the ICI status of 438 NB samples. Three ICI models were selected and ICI scores were acquired. Subgroups with high ICI scores based on immune activation signaling pathways have better overall survival. Results: The genes of immunosuppressive glycosaminoglycan biosynthesis heparan sulfate signaling pathway were markedly enriched in the low ICI score subgroup. It was inferred that compared with low ICI NB subtypes, patients with high ICI NB subtypes were more likely to respond to immunotherapy and a better prognosis. Conclusion: Notably, our ICI scores not only provided new clinical and theoretical basis for mining NB prognostic markers related to the microenvironment, but also aided new ideas for the development of new NB precision immunotherapy methods.


Introduction
Neuroblastoma (NB) is an extracranial malignant solid tumor derived from the embryonic neural crest cells, its mortality accounts for nearly 15% of all pediatric tumors [1,2]. The age of onset of NB is usually early, and the median age at diagnosis is 17 months. It usually occurs in the adrenal medulla or paravertebral ganglia, so it can locate in the abdominal cavity, thoracic cavity, neck, etc. About half of NB patients have stage distant metastases at the rst diagnosis [3]. Although the ve-year survival rate of low-risk NB is greater than 95%, the ve-year survival rate of high-risk is still less than 50% after comprehensive treatment such as surgery, radiotherapy, chemotherapy and biological therapy, because of its later recurrence, distant metastasis, and drug resistance [4].
Over that last decades, various immunotherapies have been applied to treat more and more refractory malignant cancers [5]. Immunotherapy methods related to neuroblastoma treatment have been expected to become the most promising treatment strategy, providing a new direction for its treatment. However, immunotherapeutic strategies bene t only a subset of patients, while most patients face the problem of resistance to immunotherapy [6]. With the insights into of the immune system and tumor microenvironment, new targets may facilitate NB patients for immunotherapy treatment.
High risk NB has a high degree of malignancy and is resistant to radiotherapy and chemotherapy. Therefore, new treatment strategies are needed to improve their biological behavior. The tumor-speci c antigen and immune microenvironment phenotype may be used as predictive indicators for evaluating the effect of immunotherapy. Tumor microenvironment (TME) is a dynamic system formed by the interaction of cancer cells with tumor interstitial cells, immune cells, extracellular matrix, surrounding vascular tissues, and signaling molecules [7]. In addition to the invasion, migration, and proliferation of tumor cells, tumor metastasis is also closely related to its microenvironment. The inhibitory changes and heterogeneity of TME are important factors that promote tumor progression and affect immune e cacy [8][9]. Due to multiple factors such as immunosuppressive microenvironment and lack of tumor-speci c antigen expression, the high-risk group NB is not sensitive to immunotherapy and belongs to cold tumor, while low-risk group NB has the features of hot tumor. Increased expression of tumor-associated macrophages M2 (CD163+) could decrease function of NK cells, DCs and T cells in these cold tumors [10]. The main components of in ltrating stromal cells are macrophages and broblasts. The cell-to-cell interaction between tumor-associated macrophages M2 and Cancer-associated broblasts promotes the progression of neuroblastoma by mutual induction of recruitment and activation. Studies have con rmed that these cells promote tumor-promoting activity by activating the ERK1/2 and STAT3 signaling pathways in neuroblastoma cells [11]. These cells, called CAFs-MSCs, promote the proliferation of NB cells and secrete a variety of in ammatory cytokines and chemokines to enhance resistance to chemotherapy, such as IL-6, IL-8, VEGF-A, CCL-2 And CXCL12. The Disialoganglioside (GD2) was highly expressed in NB cells, resulting in the immune function of T cells impaired and immunotherapy resistance [12]. Programmed death 1(PD-1), the main immune checkpoint receptor, expressed on human T cells, dendritic cells and natural killer T cells, while programmed death ligand 1 (PD-L1) was expressed on the surface of NB cells. The PD-1/PD-L1 pathway inhibits the proliferation and differentiation of T cells after PD-1 binds to the PD-L1 of NB cells, then induce activated T cells to transform into non-effect T cells or apoptosis [13][14].
Next generation sequencing (NGS) has become an important part of precise tumor diagnosis and treatment that used in multiple aspects such as the detection of driver genes related to tumor targeted therapy, the analysis of drug resistance mechanisms, the evaluation of cancer metastasis and prognosis, and the full dynamic monitoring of tumor patients. CIBERSORT is a tool that uses a deconvolution algorithm to predict the proportion of cells. It can analyze the cellular composition of immune tissues based on standardized expression data and energize the abundance of speci c cell types [15][16]. In this study, we explored the gene expression of TARGET-NB and GSE85047 to get a comprehensive overview of the tumor immunity by two calculation algorithms (CIBERSORT and ESTIMATE) [17]. Then we divided NB into different discrete subtypes by immune cell in ltration mode [18]. Besides, we divided NB into different discrete subtypes by immune cell in ltration mode. Establishing the ICI score to characterize various immune characteristics could promote the best choice for NB patients who respond to immunotherapy.

Neuroblastoma Data collection
The RNA data (level 3) in TARGET Data were downloaded from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET)-NB database (https://ocg.cancer.gov/programs/target/datamatrix). The Affymetrix microarray data sets GSE85047 were obtained from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) database. Clinical data such as age, INSS, survival and outcome were also downloaded from TARGET and GEO, respectively. Fragments per kilobase of transcripts per million (FPKM) data from the TARGET samples were converted to transcripts per million (TPM). In order to obtain key information from these expression matrices, Perl scripts and R software (R-project.org) were used to merge and preprocess the original data. If the gene symbol corresponds to multiple probe IDs, the average value of the probes was calculated as the expression level of the gene.
Immune score determination for the NB microenvironment and immune cellular fraction estimates We calculated the immune score of NB with the ESTIMATE algorithm. All patients were classi ed into high group and low group according to their immune/matrix score. The CIBERSORT algorithm was used to evaluate tumor in ltrating immune cell (TIIC) data in NB samples from the TARGET and GEO cohorts [19][20]. The visualization of the results was completed by using the R packages. The hierarchical agglomerative clustering of NB was performed according to the ICI mode of all patients. This analysis was executed using the R package and repeated 1000 times to ensure the stability of the classi cation.

Differentially expressed genes (DEGs)-related with the ICI clusters
Patients were classi ed into different ICI clusters base on immune-cell in ltration. Absolute fold change > 1 and p < 0.05 (after adjustment) were determined by setting the cut-off criterion for DEGs. Data was applied to analysis by package edge R [18].

Pathway and function enrichment analysis
Functional enrichment analysis of immune-cell in ltration-related DEGs was performed with the "clusterPro ler" R package to identify the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). P value < 0.05 was considered as the critical cutoff [21].

Principal component analysis (PCA)
PCA is a statistical method that applies the idea of dimensionality reduction to transform multiple variables into a few uncorrelated comprehensive variables [22]. We performed PCA to de ne the differences in in ltrating immune cells among distinct groups. The diagram of PCA could classify different in ltrating immune cells as variables to describe the differences in NB samples.

Statistical Analyzes
All statistical analyses were performed using SPSS 22.0 software. Kaplan-Meier analysis was conducted to detect the survival curve of the subgroups in the data set. The log-rank test assessed statistically

Results
The pro le of Immuno-cell In ltration in the TME of NB The study design was shown in Fig. 1. In our study, we rstly adopted the two algorithms (ESTIMATE and CIBERSORT) to analyze the immune cells subtype with R platform in NB samples. Subsequently, the ConsusClusterPlus package was performed unsupervised clustering of 438 NB samples from GSE85047 and [TARGET]-NB) to de ne different subgroups.
The Kaplan-Meier curves indicated that ICI A and ICI B were not signi cantly associated with improved outcome (log rank test, p = 0.179; Fig. 2A-C). Then we calculated the immune cell subtype in R platform between the two ICI Subtypes.
The relative content of resting dendritic cells in all samples were nearly 0. Between these two main immune subtypes, the immune cells found in the ICI cluster A mainly comprised high ImmuneScore, naive B cells, T cells regulatory (Tregs), naive CD4 T cells and CD8 T cells, while the samples in the ICI cluster B were represented with a signi cantly higher density of StromalScore, T cells CD4 memory resting, M0, M2 macrophages and Mast cells resting (Fig. 2D). In addition, the heat map of correlation coe cients of NB TME was analyzed (Fig. 2E). We found that the relative content of M2 cells in NB samples was signi cantly negatively correlated with those cells (CD8 T cells, naive B cells, T cells regulatory).

Identi ed Immune Gene Subtype
After the gene datasets from TARGET-NBL and GSE85047 samples were integrated with R (version 3.6.1), we identi ed an aggregate of 117 DEGs by limma package. Then we divided these samples into three genomes based on DEGs and called them gene clusters A-C (Fig. 3A). There were 24 DEGs signi cantly positively related to the gene cluster were named ICI gene features A, while the remaining DEGs were called ICI gene features B. Dimension reduction was used to reduce interference by Boruta algorithm. The heatmap of the 117 DEGs in R software was drawn to show the correlation between the genomic traits and clinical features (Fig. 3B) [23]. The enrichment analyses by the GO term and KEGG pathway were showed in Figs. 3C and 3D. The ICI features genes A included neutrophil activation involved in immune response, colin − 1−rich granule and neutrophil degranulation; while ICI features genes B gene were associated with extracellular matrix structural constituent, collagen-containing extracellular matrix and extracellular structure organization.
In addition, the Survival and survminer package of R software were performed to explore the relationship between survival feature and ICI gene clusters. The Kaplan-Meier analysis revealed that the samples in gene cluster A were signi cantly associated with improved outcome, while samples in remaining gene clusters had a poor prognosis (log rank test, p = 0.040; Fig. 3E). The gene cluster C was characterized by more abundant stromal cell in ltration and immunosuppression-related M2 macrophages, suggesting signi cant differences. The Wilcoxon test was performed to compare the two groups, while the Kruskal-Wallis test was applied to compare more than two groups. that the result of immunotherapy was unfavorable [24][25]. The expression of PD-L1 and CTLA4 were signi cantly different in these three genome clusters (Kruskal-Wallis, p < 0.001; Fig. 2G and 2H).
Compared with ICI gene clusters A and C, CTLA4 and PD-L1 had the lowest expression levels in ICI gene cluster B.

Construction of the ICI Score
We performed principal component analysis (PCA) to calculate the quantitative indicators of different ICI subtypes, which were labeled as ICI scores A and B, respectively. Then, we obtained prognostic-related score de ned as ICI scores and divided the patients in the TARGET and GSE85047 cohorts into high or low ICI scores subgroup. The distribution of different patients in diverse gene clusters was summary by the Alluvial graph (Fig. 4A). The Gene ICI cluster A that expressed a higher ICI scores was signi cantly related to an improved survival outcomes of 72.3% (136/188), while the samples in the gene ICI cluster B expressed a less high ICI scores witnessed a lower overall survival outcome (58.2%, 131/225). Then we analyzed the relationship between the immune tolerance of the NB patients and the prognostic value of the ICI score. A few representative immune activity-related targets (such as PD-L1, CTLA4, LAG3, IDO1 and HAVCR2) and immune checkpoint-related targets were selected [26][27]. In our study, Exception of LAG3, all selected immune activity-related targets and immune checkpoint-related targets were remarkably overexpressed in the high ICI group. The gene set enrichment analysis (GSEA) showed that the glycosaminoglycan biosynthesis heparan sulfate signaling pathway was associated with the low ICI score group, while apoptosis, cytokine receptor interaction, toll like receptor and leukocyte transendothelial migration signaling pathways were associated with the high ICI group (Figs. 4C and 4D). The Kaplan-Meier analysis revealed that the samples in high ICI score subtype were signi cantly associated with better outcome, while samples in low ICI score subtype had a unfavorable OS (p = 0.025; Fig. 4E).

Discussion
Immunotherapy have shown encouraging results and may improve the survival status and quality of high-risk NB patients [13][14]. However, the inconsistency of their exploratory study of tumor surface markers highlights the necessity of determining the ideal subgroup of NB immunotherapy, which is still a major challenge in immunotherapy. A monoclonal antibody named dinutuximab attached to GD2 on the surface of neuroblastoma cells has been performed as an immunotherapy drug for the neuroblastoma [6]. Recently, various evidence has con rmed the potential of PD-1/PD-L1 in immunotherapy of neuroblastoma [28][29]. However, only a small percentage of patients bene t from immunotherapy. Therefore, it is important to identify those subgroups that can bene t from immunotherapy. In our research, we validated a method that can quantify the integrated tumor microenvironment in NB. Our results indicated that the ICI score was not only an effective prognostic feature, but also been used to identify NB patients that could bene t from immunotherapy.
As in previous studies, patients with speci c active immunity have a signi cantly better prognosis, while high-risk NB patients with immune cell dysfunction increased immune resistance, resulting in tumor progression and unfavorable prognosis [30][31]. At rst, we analyzed the ICI from 438 NB sample cohorts in two databases and divided NB into two different immune subtypes. The results displayed that the relative content of resting dendritic cells in all samples was nearly zero. ICI cluster A was marked by high ImmuneScore, naive B cells, T cells regulatory (Tregs), naive CD4 T cells and CD8 T cells. The samples in the ICI cluster B were represented with a signi cantly higher expression of StromalScore, T cells CD4 memory resting, M0 and M2 macrophages, and Mast cells resting. Although the Kaplan-Meier survival curves was not statistically signi cant.
Studies have shown that a number of cytokines, TME components and the host's immune response have a positive impact on the anti-tumor response and maintain a dynamic balance. These molecular changes during tumorigenesis may interfere with the function of in ltrating immune cells, resulting in the destruction of the dynamic balance between immune tolerance and activity [32].
We utilized the combination of ICI and immune associated gene expression as a new method for patientspeci c customized treatment strategies. Then, the novelty immune-related genes based on the signi cant ICI gene clusters were used to explore its prognostic signi cance by integrating ICI gene clusters with survival information. The Kaplan-Meier analysis revealed that the samples in gene cluster A were signi cantly associated with improved outcome, while samples in remaining gene clusters had a poor prognosis (log rank test, p = 0.040; Fig. 2E). Patients the lowest expression of immune score and matrix score are grouped in ICI gene cluster B subtype, which are related to the immune cold phenotype. While patients with ICI gene cluster A and C subtypes showed higher in ammatory cell density and immune scores. Moreover, we found that the increased expression of M2 macrophage in ltration in ICI gene cluster C subtype was related to high stromal score, meaning an unfavorable prognosis [25].
Compared with other subtypes, samples ICI gene cluster A subtype are associated with a good prognosis, which could bene t from immunotherapy.
There is an urgent need to quantify the ICI pattern of neuroblastoma because of the heterogeneity of the individual microenvironment. Similar individual models have been well established in colorectal cancer, head and neck squamous cell carcinoma and breast cancer to improve prediction accuracy [18, [33][34].
The potential "signature subtype" was established to quantify the ICI model in our study. The Kaplan-Meier analysis showed that the samples high ICI score subtype was signi cantly associated with better outcome, while samples in low ICI score subtype had an unfavorable OS (log rank test, p = 0.025). The most immune checkpoint associated targets and immune activity associated targets were overexpressed in the high ICI subgroup. It is inferred that compared with ICI NB subtypes, patients with high ICI NB subtypes were more likely to respond to immunotherapy. For example, PD-1/PD-L1 targeted immunotherapy was more likely to be effective because of the PD-L1 expression was more expressed in the high ICI NB subgroup. The genes of immunosuppressive glycosaminoglycan biosynthesis heparan sulfate signaling pathway were remarkably enriched in the low ICI score subgroup by GSEA analysis. Destroying the integrity of the extracellular matrix and basement membrane was the primary condition for the invasion and metastasis of malignant tumors. Microscopic damage on the extracellular matrix in the tumor microenvironment is a sign of tumor aggressive growth [35]. Heparanase (HPSE), as an endogenous endoglycosidase, was highly expressed in high-risk neuroblastoma. Tumor cells secreted a large amount of heparanase to destroy the strong network structure of extracellular matrix and basement membrane, promoting the invasion and metastasis of in ammatory cells and tumor cells and related diseases Progress [36]. In the high ICI group, Apoptosis, Cytokine receptor interaction, Toll like receptor and Leukocyte transendothelial migration signaling pathways were enriched.

Conclusion
In summary, it is urgent to constantly explore novel prognostic markers for the treatment of tumor, we found that ICI score presented here is applicable in NB and can also be used for development of speci c prognostic biomarkers. We found that differences in ICI subtypes are related to individual heterogeneity and treatment e cacy. In addition, it can help develop new ideas for NB immunotherapy methods.