Jiaxing Yang1,2,3,4†
Lei Han2,3,4,5†
Yongliang Sha1,2,3,4
Yan Jin1,2,3,4Zhongyuan Li1,2,3,4
Baocheng Gong1,2,3,4Jie Li1,2,3,4Yun Liu1,2,3,4
Yangyang Wang1,2,3,4
Qiang Zhao1,2,3,4*- 1Department of Pediatric Oncology, Tianjin Medical University Cancer Institute and Hospital, Tianjin, China
- 2National Clinical Research Center for Cancer, Tianjin Medical University Cancer Institute and Hospital, Tianjin, China
- 3Key Laboratory of Cancer Prevention and Therapy, Tianjin, China
- 4Tianjin’s Clinical Research Center for Cancer, Tianjin, China
- 5Department of Cancer Molecular Diagnostics Core, Tianjin Medical University Cancer Institute and Hospital, Tianjin, China
Introduction: Gangliosides play an essential role in cancer development and progression. However, the involvement of gangliosides in the prognosis and tumor microenvironment (TME) of neuroblastoma is not entirely understood.
Methods: Consensus clustering analysis was performed to identify ganglioside-mediated molecular subtypes. LASSO-Cox analysis was conducted to identify independent prognostic genes, and a novel risk signature was constructed. The risk signature was validated internally and externally. We further explored the independent prognosis value, immune landscape, drug susceptibility, and tumor dedifferentiation of the risk signature. The role of the signature gene B3GALT4 in neuroblastoma was explored in vitro.
Results: Seventeen ganglioside-related genes were differentially expressed between INSS stage 4 and other stages, and two ganglioside-related clusters with distinct prognoses were identified. A novel risk signature integrating ten ganglioside-related prognostic genes was established. Across the train set and external validation sets, the risk signature presented high predictive accuracy and discrimination. The risk signature was an independent prognostic factor and constructed a nomogram combining multiple clinical characteristics. In the high-score group, the deficiency in antigen processing and presenting machinery, lack of immune cell infiltration, and escaping NK cells contributed substantially to immune escape. The low-score group was more responsive to immune checkpoint blockade therapy, while the high-score group showed substantial sensitivity to multiple chemotherapeutic drugs. Besides, the risk score was significantly positively correlated with the stemness index and reduced considerably in all-trans retinoic acid-treated neuroblastoma cell lines, indicating high dedifferentiation in the high-score group. Additionally, neuroblastoma cells with downregulation of B3GALT4 present with increased proliferation, invasion, and metastasis abilities in vitro.
Conclusion: The novel ganglioside-related risk signature highlights the role of ganglioside in neuroblastoma prognosis and immune landscape and helps optimize chemotherapy and immunotherapy for neuroblastoma.
1. Introduction
Neuroblastoma is the most common extracranial solid tumor in children, accounting for only 6-10% of all pediatric malignancies but 12-15% of pediatric cancer-related deaths (1–3). The prognostic heterogeneity of neuroblastoma has been wildly characterized. The 5-year event-free survival (EFS) rate in the low-intermediate risk neuroblastoma exceeds 80%, while the high-risk group, which accounts for half of the total cases, has a 5-year EFS of only 50% (4). Further survival improvement needs more precise prognostic information on neuroblastoma, and novel genetic and molecular predictive biomarkers are urgently required besides known clinical risk factors.
Gangliosides are glycosphingolipids prevalent on the surface of cells, characterized by one or more sialic acid residues on carbohydrate moieties. It is particularly in specialized membrane domains known as lipid rafts and has a role in cell adhesion and signal transduction (5–7). Gangliosides are implicated in cancer development and progression, including tumor proliferation, invasion, angiogenesis, and metastasis (8). Ganglioside GM3 could decrease the phosphorylation of epidermal growth factor receptors and inhibit the proliferation of bladder cancer (9). Ganglioside GM2 is highly expressed in pancreatic ductal adenocarcinoma and correlated with the activation of TGF-β1 signaling and the promotion of tumor invasion (10). Besides, ganglioside GM3 and GD3 are involved in angiogenesis regulation and metastasis in solid tumors (11, 12). However, the role of gangliosides in neuroblastoma is not entirely understood, and the double-edged sword function of gangliosides in regulating malignant characteristics in neuroectodermal-derived malignancies is a critical trait (8, 13, 14). The monosialogangliosides GM3 and GM1 suppress neuroblastoma, glioma, and astrocytoma proliferation by interacting with different growth factor receptors (13). By contrast, GD3 and GD2 of the b-series gangliosides contribute predominantly to tumor-promoting activities in malignancies arising from neuroectodermal cells (13). GD3 is involved in maintaining and enhancing neural stem cell and glioblastoma self-renewal abilities through EGFR activation (15, 16). Besides, GD3 and GD2 promote proliferation, motility, and invasion in various malignancies, including breast cancer, small cell lung cancer, melanoma, and osteosarcoma (13, 17–19). Importantly, GD2 has become one of the most critical tumor markers and immunotherapeutic targets for neuroblastoma. The anti-GD2 monoclonal antibody immunotherapy has been wildly conducted in neuroblastoma clinical management and presents a considerable improvement in high-risk neuroblastoma prognosis. Given the critical and complicated effects of gangliosides on tumors, it is vital to conduct in-depth studies on the role of gangliosides in neuroblastoma to predict prognosis and inform clinical management.
The remarkable efficacy of anti-GD2 monoclonal antibody suggests the tremendous potential of immunotherapy in neuroblastoma. Nonetheless, immunotherapeutic approaches in neuroblastoma continue to face several obstacles. Neuroblastoma has low immunogenicity due to low mutational load and MHC-I expression, resulting in a lack of lymphocyte infiltration and immunological activity in the tumor microenvironment (TME) (20). Additionally, various immune evasion strategies in TME can obstruct lymphocyte infiltration and activation (20). Thus, a comprehensive understanding of the TME is critical for precisely targeting neuroblastoma with immunotherapy. Interestingly, previous studies have shown that gangliosides are involved in regulating TME. Gangliosides on the surface of tumor cells or shed from cells can suppress cytotoxic T cells or dendritic cells, contributing to tumor immune evasion (7, 20). Besides, gangliosides and IFN-γ could synergistically inhibit dendritic cell activity, promoting immune suppression in the TME (21). Therefore, elucidating the role of gangliosides in TME in neuroblastoma could facilitate the understanding of tumor progression and the optimization of immune therapies.
While the involvement of gangliosides in neuroblastoma prognosis and TME remains to be explored, the bioinformatic analysis provides us with a new direction. In the present study, samples in the GSE49710 dataset were clustered based on ganglioside-related gene expression, and a risk signature was constructed to predict neuroblastoma prognosis. Additionally, we investigated the immune landscapes and escape strategies in ganglioside-related risk groups. Immunotherapy response and chemotherapeutic drug sensitivity were further explored in the high-score and low-score groups. Our findings constructed an accurate and effective prognostic signature for neuroblastoma and may help inform the treatment strategy for neuroblastoma.
2. Materials and methods
2.1. Data acquisition and preprocessing
The workflow of this study was presented in Figure S1. Expression data and corresponding clinical information were obtained from the Gene Expression Omnibus (GEO) GSE49710 (n = 498) (22) and ArrayExpress E-MTAB-8248 (n = 223) (23). The Therapeutically Applicable Research to Generate Effective Treatments (TARGET) neuroblastoma gene-expression profile and clinical data (n = 150) were acquired from the UCSC Xena database (http://xena.ucsc.edu/). The GSE49710 cohort was used to construct the risk signature, with the E-MTAB-8248 and TARGET datasets serving as external validation. The clinical baseline characteristics of three data sets were summarized in Table 1. Expression data were normalized, and log2 transformed. The expression profile and corresponding immunotherapy information for the GSE78220 cohort were retrieved from the GEO database. The expression data of neuroblastoma cell lines treated with Dimethyl sulfoxide (DMSO) or all-trans retinoic acid (ATRA) was obtained from GSE155000 in the GEO database. Thirty-four genes associated with gangliosides were identified by the Molecular Signatures Database (MSigDB; https://www.gsea-msigdb.org/gsea/msigdb) and previously published literature (Table S1) (14).
Table 1 Clinical baseline characteristics of three datasets in present study.
2.2. Consensus clustering analysis based ganglioside-related genes
A total of thirty-three ganglioside-related genes were extracted from the GSE49710 dataset. Given the significant difference between International Neuroblastoma Staging System (INSS) stage 4 and other stages (24), seventeen ganglioside-related genes were identified as differentially expressed genes (DEGs) between these two groups. Unsupervised consensus clustering analysis based on ganglioside-related DEGs was performed on GSE49710 and TARGET datasets using the “Consensus Cluster Plus” R package and the “k-means” method to discern diverse expression patterns, with the repetition number set to 1000 to ensure stability (25). The t-distributed stochastic neighbor embedding (t-SNE) analysis was conducted using the “Rtsne” package to investigate the distribution between distinct clusters.
2.3. Identification of differentially expressed genes between clusters and functional enrichment analysis
The analysis of differentially expressed genes between clusters was conducted using the “limma” R package according to the specified criteria (|log2FC| ≥1 and adjusted p-value < 0.05). The clusterProfiler R package was used to perform functional enrichment analysis on these DEGs, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses (26).
2.4. Construction and validation of the ganglioside-related risk signature
Univariate Cox regression analysis using the “survival” R package was performed to screen for potential prognostic genes within ganglioside-related DEGs in the GSE49710 dataset. The least absolute shrinkage and selection operator (LASSO)-penalized Cox regression analysis was then conducted to identify potential signature genes. Finally, the regression coefficients of ten signature genes were determined using multivariate Cox regression. The risk score was calculated by multiplying the expression value of each signature gene by its corresponding regression coefficient.
According to the median risk score, samples were classified as low-score or high-score groups. The scatter dot plots were generated to visualize the association between risk score and survival status. Survival analysis of overall survival (OS) probability was performed to evaluate the risk signature prognostic value using “survival” and “survminer” R packages. Receiver operating characteristic (ROC) curve analysis was performed by the “timeROC” R package to assess the specificity and sensitivity of the risk signature. The unique capability of risk signature was evaluated using principal component analysis (PCA) and the R package “ggplot2.” The ganglioside-related risk signature validation was carried out in GSE49710, TARGET, E-MTAB-8248, and Tianjin cohorts.
2.5. Independent prognostic analysis and construction of a nomogram
The “survival” R package was used to conduct univariate and multivariate Cox regression analysis on datasets to determine the risk signature’s predictive significance in the context of recognized prognostic indicators. Furthermore, a nomogram was constructed comprising the risk signature and several established prognostic factors by the “rms” R package to predicate 3-year and 5-year OS. Calibration plots were drawn to assess the accuracy of the nomogram.
2.6. Gene set enrichment analysis
Gene Set Enrichment Analysis (GSEA) was conducted to investigate the biological function differences in low-score and high-score groups using the R packages “clusterProfiler” and “enrichplot.” The reference gene set (“c2.cp.kegg.v7.5.1.symbols.gmt”) was acquired from MSigDB (http://www.gsea-msigdb.org/gsea/downloads.jsp).
2.7. Immune landscape of the risk signature
The low immunogenicity of neuroblastoma could be attributed to a low mutational burden and impairment in the antigen processing and presenting machinery (APM) (20). APM scores (APS) were derived in the previous study by gene set variation analysis (GSVA) based on 18 APM-related genes to represent antigen processing and presentation efficiency (Table S2) (27). As an integral element of APM, MHC-I mediates the recognition and lysis of neuroblastoma cells by cytotoxic T lymphocytes (CTL). To estimate the quantity of MHC-I activity, we obtained the gene set associated with the MHC-I protein complex pathway from the MSigDB (Table S3). Subsequently, the enrichment score (ES) was calculated in the single-sample gene set enrichment analysis (ssGSEA) using the “gsva” R package to reflect MHC-I activity, which was defined as the MHC score in this study.
The immune infiltration landscape was investigated by calculating infiltrating scores of 30 different types of TME cells by ssGSEA (Table S4) (28, 29). Besides, immune, stromal, and ESTIMATE scores were calculated by the “ESTIMATE” algorithm to reveal the distinct immune microenvironments between risk groups (30).
Immune function-related gene sets, including interferon receptor and natural killer cell cytotoxicity, were gathered from the Immunology Database and Analysis Portal (ImmPort, http://www.immport.org) database (Table S5) (31). Immune function activities were measured by ssGSEA and compared between risk groups.
2.8. Development of the stemness index
The messenger ribonucleic acid stemness index (mRNAsi) was developed using the one-class logistic regression machine learning algorithm (OCLR) based on pluripotent stem cell samples from the Progenitor Cell Biology Consortium dataset (https://www.synapse.org/, accessed on 16 January 2022). The mRNAsi had been widely used for tumor dedifferentiation and stemness prediction (32–34). The workflow was available on https://bioinformaticsfmrp.github.io/PanCanStem_Web/. The mRNAsi value was normalized to 0-1, with increased mRNAsi indicating a greater degree of dedifferentiation. The stemness index model was constructed in this work, and the mRNAsi was estimated in the GSE49710 dataset. Considering the effect of tumor purity on mRNAsi, we corrected for mRNAsi using tumor purity generated by the ESTIMATE algorithm, and the corrected mRNAsi (c_mRNAsi) was calculated as mRNAsi/tumor purity.
2.9. Drug sensitivity analysis
The half-maximal inhibitory concentration (IC50) of commonly used chemotherapeutic agents was predicted by the “pRRophetic” R package to characterize chemosensitivity in high-score and low-score groups (35).
2.10. Immunohistochemistry
A total of forty-six paraffin-embedded neuroblastoma specimens were collected at Tianjin Medical University Cancer Institute and Hospital. This study complied with the Declaration of Helsinki and was approved by the Ethics Committee of Tianjin Medical University Cancer Institute and Hospital (E20210027). Sections were deparaffinized with xylene for 30 mins and gradient concentrations of alcohol followed by rehydration. The heat-induced epitope retrieval was conducted by the Tris/EDTA buffer (Solarbio, Beijing, China), pH of 9.0, at 120°C for 3 mins, and the sections were immersed in 3% hydrogen peroxide for 30 min and incubated with the primary ST8SIA2 antibody (dilution 1:100; Rabbit polyclonal, 19736-1-AP; Proteintech), B3GALT4 antibody (dilution 1:100; Rabbit monoclonal, ab169759; abcam), and CD8 antibody (dilution 1:4000; Mouse monoclonal, 66868-1-Ig; Proteintech) at 4°C overnight. After washing with PBS and incubation with the secondary antibody (PV-6001; ZSGB-BIO; Beijing, China) at 37°C for 1 hour, the antigens were detected using DAB chromogen and counterstained with hematoxylin for 1 min. The immunoreactivity score (IRS) was generated for semi-quantitative expression and scored by two independent, experienced pathologists blinded to the clinical information. The inconsistencies were discussed to reach a unified result. The IRS considered staining intensity and the percentage of positive tumor cells. The staining intensity was assessed in four grades, including negative staining (0 points), weak staining (1 point), moderate staining (2 points), and strong staining (3 points). The percentage of positive tumor cells in the section was divided into five grades, including 0-5% (0 point), 6-25% (1 point), 26-50% (2 points), 51-75% (3 points) and 76-100% (4 points). The sample IRS was calculated by multiplying scores of the staining intensity and percentage of positive tumor cells. The precents of CD8+ T cells were quantified as the proportion of CD8A-positive cells in all cells on 200× photographs.
2.11. Cell lines and cell culture
Neuroblastoma cell lines 9464D and 975A2 were gifted from Dr. Rimas Orentas at Seattle Children’s Research Institute. Neuroblastoma cells were maintained in the high-glucose DMEM medium (Gibco) containing 10% FBS (BI) and 1% penicillin/streptomycin (Gibco). The cells were cultured in a humidified incubator at 37°C in a 5% CO2 atmosphere.
2.12. Cell transfection
The small interfering RNA (siRNA) targeting B3GALT4, and negative control siRNA (si-NC) were purchased from General Biol (Anhui, China). The manufacturer’s instruction was followed for cell transfection. Neuroblastoma cells were transfected with siRNA using the transfection reagent Lipofectamine®2000 (Invitrogen). The transfected cells were collected for further experiments after 24h. Quantitative realtime PCR assay and western blot analysis were performed to verify the knockdown efficiency.
2.13. Western blot analysis
Cells were lysed in RIPA lysis buffer (Solarbio) for protein extraction, and the protein concentrations were evaluated by the BCA method. The proteins were separated by 10% SDS-polyacrylamide gel electrophoresis (PAGE) and transferred to PVDF membranes. After incubation in 5% skimmed milk for 1 h at room temperature, the membranes were incubated overnight at 4°C with primary antibodies against B3GALT4 (dilution 1:1000; Rabbit monoclonal, ab169759; abcam). After incubation with the secondary antibody and wash with TBS-T three times, the band images were visualized by the enhanced chemiluminescence kit.
2.14. Real-time quantitative PCR
The total RNA was extracted from cells by Trizol reagent (Invitrogen) and converted to cDNA using the PrimScript RT Master Mix (Takara). cDNA amplification was carried out by SYBR Green PCR Kit (Takara) according to the program: 5 seconds at 95°C for the denaturation, 34 seconds at 60°C for annealing, followed by 30 seconds at 72°C for extension, and forty cycles were completed in total. The primer sequences were designed as follows: B3GALT4: F: 5’-AACGCCATTCGGGCATCTT-3’, R: 5’-GTTGCGGTAGGAATCCTGGAA-3’; GAPDH: F: 5’-ACCCTTAAGAGGGATGCTGC-3’, R: 5’-CCCAATACGGCCAAATCCGT-3’. The 2−ΔΔCt value was employed to quantify the relative gene expression levels with GAPDH as the endogenous control.
2.15. Cell proliferation and colony forming assay
Cell viability was measured by the Cell Counting Kit-8 (CCK8) assay. Neuroblastoma cells were plated into 96-well plates for 24h, 48h, and 72h, followed by adding 100µ CCK8 solution (Solarbio) and incubating for 2 hours at 37°C. The absorbance of each well was measured at a wavelength of 450 nm (OD450) with a microplate reader. Cells were planted and cultured for 2 weeks in each well of a 6-well plate for cell colony formation assay. The colonies were fixed with 4% paraformaldehyde for 15 minutes, stained with 0.1% crystal violet for 20 minutes at room temperature, and quantified by ImageJ software after being photographed.
2.16. Cell invasion and migration assay
The transwell assay was performed to evaluate cell invasion and migration ability. Eight-micrometer pore-size transwell filters (Corning) were put in a 24-well plate for the migration assay, while the upper chambers plated with matrigel (BD Biosciences) for the invasion assay. Cells in 200 μl FBS-free medium were seeded onto each upper chamber, and the lower chamber was added with 600 μl medium with 10% FBS. After being cultured for 24 hours at 37°C, these invasive and metastatic cells in the lower side of the filter were fixed by 4% paraformaldehyde, stained with 0.1% crystal violet solution, and photographed.
2.17. Sample collection for RNA sequencing in the Tianjin cohort
Twenty-six neuroblastoma biopsies and corresponding clinical information were collected at Tianjin Medical University Cancer Institute and Hospital. This study complied with the Declaration of Helsinki and was approved by the Ethics Committee of Tianjin Medical University Cancer Institute and Hospital (E20210027).
2.18. RNA quantification and qualification
RNA quantification and qualification were performed according to the following steps. RNA purity and concentration were generated by NanoDrop 2000, and RNA integrity and quantity were quantified by the Agilent 2100/4200 system.
2.19. Library construction
The messenger RNA was extracted from total RNA and fragmented into 300-350 bp fragments. The reverse transcription was conducted using fragmented RNA and dNTPs (dATP, dTTP, dCTP, and dGTP) to synthesize the first strand cDNA, followed by the synthesis of the second strand cDNA. After the double-strand cDNA remaining overhangs were converted into blunt ends by exonuclease/polymerase, 3’ ends of DNA fragments were adenylated, and sequencing adaptors were ligated to the cDNA. Subsequently, the library fragments were purified. The PCR was used to amplify the template, and the product was purified to form the final library.
2.20. Sequencing and quality control of the raw data
After library preparation and sample pooling, Illumina sequencing was performed on the samples. Raw data in the formation of FASTQ were processed through in-house perl scripts. Clean data were formed by reads without low-quality or adapter and ploy-N. The clean data’s Q20, Q30, and GC content were assessed. The clean reads were mapped to the silva database to eliminate the rRNA.
2.21. Reads mapping and quantification of gene expression level
Paired-end clean reads were aligned to the reference genome (hg19) using Hisat2 (36). Featurecount was used to count the reads numbers mapped to each gene (37).
2.22. Statistical analysis
All statistical analysis was performed through R software (version 4.1.2). Survival curves were generated by the Kaplan-Meier method and log-rank test for statistical tests. Spearman rank correlation was used to analyze the correlations between continuous variables. The Mann-Whitney Wilcoxon or Kruskal-Wallis test was used to compare continuous variables between groups. The Pearson chi-square test was used to compare categorical variables across groups. Two-sided p < 0.05 was considered statistically significant.
3. Results
3.1. Consensus clustering analysis of ganglioside-related genes identified two clusters of neuroblastoma with different outcomes
The International Neuroblastoma Staging System (INSS) had been wildly used in neuroblastoma clinical management, and INSS stage 4 was an independent risk factor for neuroblastoma (4). Besides prognosis, there were also significant differences in biological characteristics between stage 4 and other stages. To screen crucial ganglioside-related genes in neuroblastoma, we identified seventeen ganglioside-related DEGs between the INSS stage 4 (high stage) and other stages (low stage) (Figure 1A). Consensus clustering analysis was performed to classify patients with distinct ganglioside-mediated patterns in GSE49710 based on seventeen ganglioside-related DEGs expression. k = 2 was selected as the ideal option for cluster construction, and 498 samples were allocated to clusters A and B, with 222 and 276 samples, respectively (Figure 1B, Figure S2; Table S6). As shown in Figure 1C, the two clusters could be clearly distinguished in the t-SNE analysis. Survival analysis revealed the significant survival advantage of cluster B over cluster A (P=0.005, Figure 1D). Additionally, the TARGET dataset was used to verify the consensus clustering result. Two distinct clusters with significantly different prognoses were identified in the TARGET dataset (Figure S3A-M; Table S7), suggesting the stability of the clustering result.
Figure 1 Identification of ganglioside-related clusters in the GSE49710 dataset. (A) The expression of ganglioside-related genes that were differentially expressed between samples with high stage (INSS stage 4) and low stage (other stages). (B) Identification of two ganglioside-related clusters according to the consensus clustering matrix (k = 2) in the GSE49710 cohort. (C) The t-SNE analysis revealed a clear distinction between the two clusters. (D) Kaplan–Meier curves of overall survival (OS) in the GSE49710 cohort between different clusters.
The analysis of the differentially expressed genes was performed between two clusters in the GSE49710 dataset to explore the further difference in the biological function in identified clusters. One hundred and eight DEGs were finally identified according to |logFC| > 1 and adjusted p-value < 0.05 (Table S8). Interestingly, GO functional enrichment analysis showed these DEGs were significantly enriched in neural crest cell development of biological processes, the postsynaptic membrane of cellular components, and signaling receptor activator activity of molecular function (Figure S4). In the KEGG functional enrichment analysis, DEGs were enriched in the neuroactive ligand receptor interaction and cAMP signaling pathways (Figure S4). These results suggested a potential role for gangliosides in neuroblastoma differentiation. Consistent with the functional enrichment in receptor ligand activity, the role of gangliosides in cancer cell signaling had been widely characterized (18).
3.2. Development and validation of the novel ganglioside-related risk signature
?>Given the prognostic significance of different clusters, a ganglioside-related risk signature based on two clusters was developed in the GSE49710 cohort to predict individual prognosis accurately. Firstly, the univariate Cox regression analysis of ganglioside-related DEGs was conducted, and sixteen potential prognostic genes were recognized (Figure 2A; Table S9). Then, Lasso-penalized Cox regression analysis was performed, followed by multivariate Cox regression to identify ten independent prognostic genes and corresponding regression coefficients. These genes included ABCA2, B3GALT4, NEU4, ST3GAL1, ST3GAL3, ST6GALNAC4, ST6GALNAC5, ST6GALNAC6, ST8SIA2, and ST8SIA3 (Figures 2B, C). These ten genes were used to establish the ganglioside-related risk signature, and the risk scores were derived from the expression values of each signature gene and its corresponding regression coefficient. As shown in Figure 2D, cluster A presenting with a poor prognosis, received significantly higher risk scores than cluster B. Samples were then divided into high-score and low-score groups based on the median risk score. The expression of ganglioside-related genes was substantially different between the two groups, with eight elevated genes in the high-score group (Figure 2E).
Figure 2 Construction of the ganglioside-related risk signature. (A) The forest map of ganglioside-related differentially expressed genes (DEGs) in GSE49710 cohort generated by univariate Cox analysis (P < 0.05). (B,C) LASSO Cox regression analysis of ganglioside-related prognostic DEGs. (D) The distribution of risk scores between two ganglioside-related clusters. (E) The expression of ganglioside-related genes between the low-score and high-score group. (*: P < 0.05, **: P < 0.01, ***: P < 0.001).
Neuroblastoma with high expression of the b-series gangliosides, including GD1b and GT1b, typically presented an excellent prognosis (14). Interestingly, B3GALT4 catalyzed the first step in converting GD2 to more complex b-series gangliosides. Considering B3GALT4 serves as the connecting link of ganglioside in neuroblastoma, we further performed immunohistochemistry to validate the expression of B3GALT4 in clinical neuroblastoma specimens. Consistent with the results in the GSE49710 dataset, the immunohistochemistry analysis showed a significantly low expression of B3GALT4 in samples with a high stage (Figures 3A, B). ST8SIA2 was involved in the developmental regulation of polysialic acid and modulated neuroblastoma adhesion and metastasis. Interestingly, mRNA levels of ST8SIA2 were highest in stages 1 and 4s neuroblastoma (38). Consistent with this result, our study also showed that the expression level of ST8SIA2 was significantly upregulated in low-stage samples (Figures 3C, D).
Figure 3 Expression validation of ganglioside-related genes by immunohistochemistry (IHC). (A) Representative IHC images showing the expression of B3GALT4 in different stages. Magnification, ×200, ×400. (B) The comparison of B3GALT4 immunoreactivity score (IRS) between samples with high stage and low stage. (C) Representative IHC images showing the expression of ST8SIA2 in different stages. Magnification, ×200, ×400. (D) The ST8SIA2 immunoreactivity score (IRS) comparison between samples with high stage and low stage.
Besides, these signature genes were interconnected (Figure 4A), and the corresponding regression coefficients were presented in Figure 4B and Table S10. The risk signature was validated in the GSE49710, E-MTAB-8248, and TARGET cohorts. Samples in the GSE49710 cohort were divided into two groups based on the median risk score, and an increase in risk score was associated with a decrease in survival time (Figure 4C). As expected, those samples classified as the high-score group had a considerably poorer prognosis than samples classified as the low-score group (Figure 4D). Additionally, ROC analysis demonstrated that the area under the curve (AUC) values of the risk signature for 3- and 5-year OS prediction were 0.891 and 0.902, respectively (Figure 4E). The risk signature performed better in predicting 3-year OS than established clinical prognostic factors (Figure 4F).
Figure 4 Validation of the ganglioside-related risk signature in the GSE49710 dataset. (A) The correlation coefficient between risk score and signature genes. (B) Coefficient values of ten signature genes. (C) The distribution of the risk score (top) and survival time (bottom) in the GSE49710 cohort. (D) Kaplan–Meier curves of overall survival (OS) in GSE49710 cohort between different risk groups. (E) The receiver operating characteristic (ROC) curves of the risk signature for 3-year and 5-year OS prediction in the GSE49710 cohort. (F) The ROC curves of clinical prognostic factors for 3-year OS prediction in the GSE49710 cohort.
Moreover, validation was conducted on the E-MTAB-8248 and TARGET cohorts. Consistent with the train set, the increase in risk score was accompanied by a decrease in survival time in E-MTAB-8248 (Figure S5A) and TARGET (Figure S5E) cohorts. The risk signature could accurately predict 3-year OS in E-MTAB-8248 (Figure S5B) and TARGET (Figure S5F) cohorts, with AUC values of 0.807 and 0.667, better than clinical characteristics (Figure S5C and S5G). Importantly, the risk signature had a vital prognosis predictive value in both cohorts, exhibiting significantly poor prognosis in the high-score group (Figure S5D and S5H). The PCA analysis suggested that the high-score group could also be separated from the low-score group in all cohorts (Figure S5I-K). In addition, even in children in COG high-risk group or older than 18 months, the prognosis of samples with a high score was significantly worse than that of low-score samples (Figures S6A-D).
Additionally, we validated the reliability of risk signature with tissue samples in our center. A total of twenty-six samples with RNA-seq data were included in this study, and the corresponding clinical characteristics were shown in Tables S11–12. As expected, high-score samples in the Tianjin cohort presented low survival time (Figure 5A). The risk signature performed excellent sensitivity and specificity in predicting 3-year and 5-year OS (Figure 5B). The overall survival and event-free survival were significantly worse for the sample in the high-score group compared with low-score samples (P < 0.05, Figures 5C, D). In conclusion, the risk signature could effectively predict neuroblastoma prognosis after thorough evaluation and validation.
Figure 5 Validation of the ganglioside-related risk signature in the Tianjin cohort. (A) The distribution of the risk score (top) and survival time (bottom) in the Tianjin cohort. (B) The ROC curves of the risk signature for 3-year and 5-year OS prediction in the Tianjin cohort. (C, D) Kaplan–Meier curves of overall survival (C) and event-free survival (D) in Tianjin cohort between different risk groups.
3.3. Clinical correlation analysis, independent prognosis analysis, and construction of a nomogram
There had been several established prognostic factors in neuroblastoma, including age, INSS stage, MYCN status, and the clinical risk classification system. We explored the correlation between the risk signature and these prognostic factors. As illustrated in Figure 6A, the risk score was significantly associated with age, MYCN amplification, clinical risk group, stage, and progression. Patients with unfavorable clinical characteristics present high scores (Figures 6B-F). Similar results were observed in the Tianjin cohort, and there were significant high scores in samples with advanced stage and unfavorable histology (Figure S7).
Figure 6 Correlation analysis between clinical characteristics and the risk signature in GSE49710 cohort. (A) Correlation analysis between the risk signature and clinical characteristics in GSE49710 cohort. (B-F) The comparison of risk scores between samples with different clinical characteristics, including INSS stage (B), age (C), MYCN status (D), progression (E), and COG risk groups (F). (***: P < 0.001).
In light of the correlation between the risk signature and clinical characteristics, we investigated the independent prognostic significance of the risk signature. After univariate and multivariate cox regression analysis, the risk signature was identified as an independent prognostic factor in the GSE4910 (Figures 7A, B), E-MTAB-8248 (Figure S8A), and TARGET (Figure S8B) cohorts, respectively. To optimize the clinical utilization in individual prognosis prediction, we incorporated the risk signature and several clinical risk factors to construct a nomogram in the GSE49710 cohort (Figure 7C). The nomogram could assign a score to each prognostic factor and predict 3-year and 5-year OS based on the sum of scores in each sample (Figure 7C). The calibration curves were plotted to evaluate the accuracy of the nomogram, and the nomogram prediction curves were quite close to standard curves in GSE49710 (Figure 7D), E-MTAB-8248 (Figure 7E), and TARGET (Figure S8C) cohorts, suggesting an excellent accuracy for prognosis prediction in all datasets. The nomogram incorporated the risk score and multiple established prognostic factors and could precisely predict neuroblastoma prognosis.
Figure 7 Independent prognosis analysis of the risk signature and construction of a nomogram for 3-year and 5-year overall survival (OS) prediction. (A) The univariate Cox regression analysis in GSE49710 cohort (B). The multivariate Cox regression analysis in GSE49710 cohort. (C) The establishment of a nomogram that predicted 3-year and 5-year OS in the GSE49710 cohort. (D, E) Calibration curves of the nomogram in the prediction of 3-year and 5-year OS in GSE49710 (D) and E-MTAB-8248 (E) cohorts. *: P < 0.05, **: P < 0.01, ***: P < 0.001.
3.4. Immune landscape of the ganglioside-related risk signature
The gene set enrichment analysis (GSEA) was conducted to elucidate the biological functions behind the variations in prognosis between low-score and high-score groups. As shown in Figure S9A, the high-score group was significantly enriched in tumorigenic pathways, including cell cycle, DNA replication, homologous recombination, ribosome, and spliceosome. Interestingly, the low-score group was enriched in antigen processing and presentation and cell adhesion pathways, implying underlying immune landscape differences between the two groups (Figure S9B).
Impairment of the antigen-presenting machinery (APM) contributed heavily to the low immunogenicity of neuroblastoma. It was widely accepted that antigen presentation through MHC-I molecules did not function in neuroblastoma due to low expression levels (20, 39). These characteristics made neuroblastoma cells almost undetectable to CD8 T cells. Interestingly, the low-score group was enriched in the antigen processing and presentation pathway, implying a potential role of ganglioside in APM of neuroblastoma.
We introduced the APM scores (APS) developed in previous research (detailed in “Materials and Methods”) as a proxy for antigen processing and presentation efficiency. As seen in Figure 8A, the APS was significantly and adversely linked with the risk score, showing that APM in the high-score group was suppressed (Figure 8A). Considering MHC-I molecules as a critical component of APM, the MHC score (described in “Materials and Methods”) was developed using ssGSEA to assess MHC-I activity. As expected, there was a strong negative association between risk score and MHC score as well as activated CD8 T cell infiltration (Figures 8B, C), and increased MHC-I activity was associated with high activated CD8 T cell infiltration (Figure 8D). Additionally, consistent with the concept that IFNγ could increase MHC-I expression in neuroblastoma (40), the MHC score was substantially and positively correlated with the activity of interferon receptors (Figure 8E). In general, APM impairment contributed significantly to immune escape in the high-score group, and the low activity of MHC-I may be an important cause.
Figure 8 The immune landscape of the risk signature. (A) The comparison of antigen-presenting machinery (APM) score between different risk groups (top) and correlation analysis between the risk score and the APM score (bottom). (B) The comparison of MHC score between different risk groups (top) and correlation analysis between the risk score and the MHC score (bottom). (C) The comparison of activated CD8 T cell infiltration between different risk groups (top) and correlation analysis between the risk score and the activated CD8 T cell infiltration (bottom). (D, E) The correlation analysis between the MHC score and the activated CD8 T cell infiltration (D) as well as interferon receptor activity (E). (F) Comparison of the stromal, immune, and ESTIMATE scores between high-score and low-score groups. (G) Different infiltration levels of immune infiltrating cells between low-score and high-score groups. (*:P< 0.05, **: P < 0.01, ***: P < 0.001).
Moreover, the absence of leukocytes and the presence of immunosuppressive myeloid and stromal cells were also efficient strategies for tumor immune evasion in neuroblastoma. Interestingly, ganglioside soluble GD2 could reduce T cell proliferation, suggesting the role of ganglioside in tumor immune evasion of neuroblastoma. The ESTIMATE method was used to determine the infiltration levels of immune cells and stromal cells. The high-score group presented significantly lower stromal and immune scores than the low-score group (Figure 8F). As illustrated in Figure 8G, most immune infiltrating cells were infiltrated at low levels in the high-score group, demonstrating an immune escape mechanism resulting from the lack of immune cells in the high-score group. In summary, the immune escape strategy in the high-score group was partially attributed to the lack of immune cell infiltration.
Tumor-infiltrating lymphocytes were essential constituents of the tumor immune microenvironment in neuroblastoma. Recent studies revealed reduced CD8+ T lymphocyte infiltration in high-risk and advanced-stage neuroblastoma (20). Children with a rising rate of CD8+ T lymphocytes had a better prognosis, highlighting that strengthening CD8+ T-cell responses would be a promising therapy opportunity (41). Consistent with previous studies, our results showed that samples with high CD8+ T lymphocyte infiltration presented an excellent prognosis in the GSE49710 cohort (Figure 9A). Furthermore, the risk signature was a potent indicator of CD8+ T lymphocyte infiltration. Compared to samples with a high score, samples with a low score displayed a significantly greater infiltration of CD8+ T-cells (Figure 9B). The risk score was adversely and significantly associated with CD8+ T-cell infiltration and CD8A expression (Figure 9C). We further performed immunohistochemistry on tissue samples from the Tianjin cohort to confirm this finding. As indicated in Figures 9D, E, there was a considerable decrease in CD8+ T-cell infiltration in the high-score group compared to the low-score group (Figures 9D, E). In summary, the risk signature could reliably predict CD8+ T-cell infiltration in the immune microenvironment of neuroblastoma.
Figure 9 The high-score group presented low CD8+ T-cell infiltration. (A) Kaplan–Meier curves of overall survival between samples with different infiltration of CD8+ T-cell in the GSE49710 dataset. (B) The comparison of CD8+ T-cell infiltration between high-score and low-score groups in the GSE49710 dataset. (C) The correlation analysis of the risk score, CD8+ T-cell infiltration, CD8A expression and immune-score in the GSE49710 dataset. (D) Representative immunohistochemistry (IHC) images showing the infiltration of CD8+ T-cell between different risk groups in Tianjin cohort. Magnification, ×200, ×400. (E) The comparison of CD8+ T-cell infiltration between different risk groups in Tianjin cohort.
It had been shown that gangliosides were involved in regulating NK cell cytotoxicity through multiple mechanisms (42, 43). Importantly, antibody-dependent cell-mediated cytotoxicity (ADCC) mediated by NK cells and neutrophils was a critical mechanism for anti-GD2 antibody efficacy (44). Neuroblastoma cells lacking MHC-I molecules, which served as the ligands for killer inhibitory receptors, should be particularly susceptible to NK cells. Notably, the high-score group with reduced MHC-I activity exhibited limited NK cell infiltration and cytotoxicity (Figures 10A, B). NK-mediated elimination of neuroblastoma seemed to be shielded by other mechanisms that modify the balance of activating and inhibitory signals on NK cells (20). Therefore, we compared the expression of ligands for NK-activating receptors between two risk groups. Apart from PVR, ligands for the NK cell-activating receptors DNAM-1 and NKG2D, such as MICA, MICB, and ULBP1, were downregulated in the high-score group (Figure 10C). The downregulation of NK-activating receptors may be the potential reason for the inhibition of NK cells. Taken together, escaping NK cells was an essential driver of immunosuppression in the high-score group.
Figure 10 Comparison of NK cell activity and immune checkpoints in the risk signature (A) The comparison of the NK cell infiltration level between different risk groups (top) and correlation analysis between the risk score and the NK cell infiltration level (bottom). (B) The comparison of the NK cell cytotoxicity between different risk groups (top) and correlation analysis between the risk score and the NK cell cytotoxicity (bottom). (C) The boxplot exhibited different expression levels of ligands for NK cell-activating receptors between low-score and high-score groups. (D) The correlation analysis between the expression of signature genes and established immune checkpoints. (E) The boxplot demonstrated different expression levels of immune checkpoints between low-score and high-score groups. (*: P < 0.05, **: P < 0.01, ***: P < 0.001).
Immune checkpoints (ICs) in the tumor microenvironment could significantly impact the reactivity of tumor-infiltrating lymphocytes to neuroblastoma (20). Signature genes were significantly associated with several ICs, indicating that ICs may play a role in the ganglioside-related risk signature (Figure 10D). As seen in Figure 10E, the low-score group had increased expression of PD-L1 (CD274), CD200R1, and CD200, while the high-score group had increased expression of B7-H3 (CD276). Briefly, immune evasion strategies were mediated by different ICs in low-score and high-score groups.
3.5. The ganglioside-related risk signature could predict immunotherapeutic response and chemotherapy sensitivity
Immune checkpoint inhibitors (ICIs) had revolutionized cancer treatment. However, only a small number of patients were responding (45). Considering the dramatic differences in the immune landscape between the two groups, we investigated the risk signature’s predictive ability for immunotherapeutic benefits. Due to the lack of expression data for immunotherapy of neuroblastoma, we used the immunotherapy dataset of melanoma, which is also a neuroendocrine tumor, to investigate the role of the risk signature in immunotherapeutic response prediction. Figure 11A depicted the distribution of treatment response at different risk scores in GSE78220, and the risk score was significantly lower in the immunotherapy-responsive group than in the non-responsive group (Figure 11B). Additionally, samples with high scores had a poor prognosis (Figure 11C). The risk signature had excellent discrimination in predicting immunotherapy response, presenting an AUC of 0.728 (Figure 11D). In short, our findings showed that the ganglioside-related risk signature could accurately predict responsiveness to immune checkpoint blockade therapy. The low-score group showed a better response to immunotherapy than the high-score group.
Figure 11 The ganglioside-related risk signature could predict immunotherapeutic response and chemotherapy sensitivity. (A) The distribution of immunotherapy responses at different risk scores in the GSE78220 cohort. (B) The comparison of risk scores between non-response and response groups in the GSE78220 cohort. (C) Kaplan–Meier curves of overall survival in GSE78220 cohort between different risk groups. (D) The receiver operating characteristic (ROC) curves for immunotherapy response prediction in the GSE78220 cohort. (E) The boxplot demonstrated different IC50 (the half maximal inhibitory concentration) values of four chemotherapeutic drugs, including cisplatin, doxorubicin, etoposide and vinblastine, between low-score and high-score groups. (***: P < 0.001).
Chemotherapy was the cornerstone of neuroblastoma treatment, and we examined the susceptibility of the two groups to commonly used chemotherapeutic drugs in the established treatment regimen. The IC50 values of four chemotherapeutic drugs were compared between two groups: cisplatin, doxorubicin, etoposide, vinblastine. Interestingly, the IC50 values for these chemotherapeutic drugs were significantly lower in the high-score group (Figure 11E). Function enrichment in the cell cycle and DNA replication of the high-score group may be the potential reason (Figure S9). These results suggested that the high-score group could still benefit from the established chemotherapy regime.
3.6. The high-score group presented with a high degree of dedifferentiation
The previous study showed that ganglioside was implicated in maintaining neural stem cell self-renewal capacity (15). Furthermore, the critical implications of gangliosides in tumor stem cells have been frequently highlighted (16, 46–48). Recent research suggested that the low MHC-I expression of neuroblastoma may reflect the undifferentiated state of the neural crest (20). These findings implied that gangliosides might contribute to the undifferentiated state of neuroblastoma. Interestingly, DEGs between ganglioside-related clusters were highly enriched in neural crest cell development (Figure S4). Therefore, we explored the role of ganglioside-related risk signature in reflecting the degree of neuroblastoma dedifferentiation. Firstly, we investigated the connection between the risk score and the corrected stemness index mRNAsi (c_mRNAsi). Figure 12A showed a significant positive correlation between the risk score and c_mRANsi (R = 0.62, P < 0.001). Similarly, all signature genes were significantly correlated with the c_mRNAsi (Figure 12B). We subsequently validated this result in a dataset of neuroblastoma cell lines. It had been known that all-trans retinoic acid (ATRA) could induce differentiation in both primary neuroblastomas and cell lines. As shown in Figure 12C, neuroblastoma cells treated with ATRA exhibited a reduced risk score in both BE2C and NGP cell lines. In addition, the risk score was significantly positively associated with multiple stemness markers of neuroblastoma, including CD133, EZH2, and OCT4 (Figures 12D-F). These results indicated that the high-score group presented with a high degree of dedifferentiation.
Figure 12 The high-score group presented with a high degree of dedifferentiation. (A) The comparison of the corrected mRNAsi (c_mRNAsi) between different risk groups (top) and correlation analysis between the risk score and the c_mRNAsi (bottom). (B) The correlation analysis between the c_mRNAsi and signature gene expression. (C) The comparison of risk scores between neuroblastoma cells treated with dimethyl sulfoxide (DMSO) or all-trans retinoic acid (ATRA) in BE2C and NGP cell lines. (D-F) The risk score was significantly positively associated with multiple stemness markers of neuroblastoma, including CD133 (D), EZH2 (E), and OCT4 (F). (*: P < 0.05, **: P < 0.01, ***: P < 0.001).
3.7. Downregulated B3GALT4 promoted the progress of neuroblastoma cells
Considering B3GALT4 as the connecting link of ganglioside and with the highest absolute value of regression coefficient in signature genes, the siRNA of B3GALT4 was transfected into 9464D and 975A2 cells, and the role of B3GALT4 in neuroblastoma was explored. The western blot and RT-qPCR were performed to verify the downregulation of B3GALT4 expression after transfection for further experiments (Figures 13A, B). As shown in Figure 13C, the CCK-8 assay demonstrated the knockdown of B3GALT4 significantly promoted cell proliferation in both cell lines (Figure 13C). Consistently, colony formation assay showed that neuroblastoma cells with downregulation of B3GALT4 present more cell clones than the control group in both cell lines (Figures 13D, E). Besides, cells detected in the lower chamber were significantly increased after transfected with si-B3GALT4 compared to the control group in both invasion and migration assays (Figures 13F, G). Our findings suggested that B3GALT4 could inhibit the progression of neuroblastoma, verifying the protective role of B3GALT4 in the risk signature.
Figure 13 The downregulation of B3GALT4 could promote the progression of neuroblastoma. (A, B) The western blot analysis (A) and quantitative real-time PCR (B) was performed to validate the downregulation of B3GALT4 after transfection with siRNA in 9464D and 975A2 cells. (C) The CCK-8 assay was performed to measure the proliferation capacity of 9464D and 975A2 cells. (D, E) The colony formation assay (D) and corresponding statistical analysis (E) of 9464D and 975A2 cells. (F, G) The transwell assays were conducted to determine the effect of down-regulated B3GALT4 on neuroblastoma migration and invasion capacity in 9464D (F) and 975A2 cells (G). (***, P< 0.001; ****, P< 0.0001).
4. Discussion
Neuroblastoma is an extraordinarily lethal childhood tumor characterized by high heterogeneity. Precise and efficient prognostic prediction is critical to guide treatment. In the present study, we identified two ganglioside-related clusters with differential expression patterns and outcomes, indicating a substantial prognostic significance of gangliosides in neuroblastoma. To effectively predict the individualized prognosis of neuroblastoma based on gangliosides, we conducted the Lasso-penalized Cox regression analysis on ganglioside-related DEGs to identify independent prognostic genes and developed a ten-gene risk signature. The risk signature showed excellent discrimination and accuracy in GSE49710, TARGET, E-MTAB-8248, and Tianjin cohorts. Additionally, the risk signature was significantly related to several previously identified prognostic markers and was demonstrated to be an independent prognostic factor for neuroblastoma. Besides, a nomogram incorporating multiple established clinical prognostic characteristics was developed and verified. In summary, we developed a novel ganglioside-related risk signature in neuroblastoma that enables reliable and individualized prognosis prediction.
There are ten genes in the ganglioside-related risk signature. However, few of them are identified as prognostic genes in neuroblastoma. ABCA2 encodes a membrane-associated protein belonging to the ATP-binding cassette transporter superfamily and is involved in the metabolism of gangliosides (49). Additionally, ABCA2 is overexpressed in pediatric acute lymphoblastic leukemia and may contribute to multidrug resistance (50). B3GALT4 is involved in synthesizing GM1/GD1 gangliosides and has been identified as a prognostic marker for osteosarcoma and neuroblastoma (51, 52). Our results also showed the inhibition of B3GALT4 could significantly increase neuroblastoma cell proliferation, migration, and invasion in vitro. Notably, NEU4 has been identified as a potential regulator of neuronal development, with overexpression promoting the acquisition of a stem cell-like phenotype in neuroblastoma cells (53). The sialyltransferases are required to synthesize gangliosides, and their aberrant expression is closely related to a poor prognosis in tumors (54, 55). Among them, ST3GAL1 overexpression promotes epithelial-mesenchymal transition, migration, and invasion in ovarian cancer (56), and ST3GAL3 downregulation inhibits pancreatic cancer cell migration and invasion (57). Additionally, the sialyltransferases ST6GalNAc4, ST6GalNAc5, and ST6GalNAc6 contribute to the synthesis and metabolism of the gangliosides GD1a and GM1b. Overexpression of ST6GalNAc4 has been crucial for tumor cell glycosylation modification and lung cancer metastasis, although the roles of ST6GalNAc5 and ST6GalNAc6 in malignancies remain unclear (55). Furthermore, the polysialyltransferase ST8SIA2 is also implicated in small cell lung cancer and glioma metastasis and invasion (55). Consistent with ST8SIA3 as a risk factor in the present study, it mediates the sialylation of GM3 and GD3 and promotes survival, proliferation, clonogenicity, and migration of glioblastoma cells (55). To summarize, these signature genes are intimately engaged in the synthesis and metabolism of gangliosides and play a critical role in tumor development and progress, supporting the predictive value of the risk signature in neuroblastoma.
Our result revealed that the low-score group was enriched in the antigen processing and presentation pathway, implying potential differences in the immune microenvironment between the two groups. Gangliosides have been identified as potent inhibitors of the cellular immune response. Soluble GD2 shed from neuroblastoma cells has been shown to suppress T cell proliferation and contribute to tumor immune evasion in neuroblastoma (8, 13, 20). Moreover, anti-GD2 antibodies have been introduced into the standard treatment regime for high-risk neuroblastoma. Given the critical role of gangliosides in the neuroblastoma immune microenvironment, we investigated and compared the immune landscape between the low-score and high-score groups.
Neuroblastoma has low immunogenicity, characterized by a low mutational burden and abnormalities in the antigen processing and presenting machinery. Our findings indicated that the high-score group exhibited an impairment in APM and low MHC-I activity, suggesting that the high-score group was less immunogenic than the low-score group. Interestingly, increasing researches indicate that the fundamental reason for MHC-I suppression may be the embryonic origin of neuroblastoma. Neuroblastoma seems to represent the underdeveloped neural crest state, characterized by low MHC-I expression. Neuroblastoma cell line differentiation is associated with increased MHC-I expression (20, 58). In the present study, low MHC-I expression in the high-score group may represent a substantial degree of dedifferentiation. High c_mRNAsi in the high-score group supports this hypothesis. We also found that ATRA-induced differentiated neuroblastoma cell lines exhibited reduced scores. The risk score was also significantly positively associated with multiple stemness markers of neuroblastoma. Generally, these results suggested that gangliosides may play a role in the formation and dedifferentiation of neuroblastoma, and the ganglioside-related risk signature could reflect the degree of neuroblastoma dedifferentiation.
Low immunogenicity leads to insufficient infiltration of lymphocytes into the tumor and poor anti-tumor reactivity (20). The high-score group with low immunogenicity showed significantly reduced immune scores and infiltration levels. The CD8+ T-cell, a crucial component of the immune response to tumors, has been identified as one of the most significant immunotherapy targets for tumors (59). Recent research has demonstrated that CD8+ T-cell infiltration is substantially related to neuroblastoma prognosis (20). We discovered that the risk signature could predict CD8+ T-cell infiltration accurately and validated it using tissue samples in the Tianjin cohort. Interestingly, several immunosuppressive and stromal cells were substantially infiltrated in the low-score group. ST3GAL1 has been implicated in transforming tumor-associated macrophage differentiation to a more suppressive phenotype (55). Consistently, our results also showed that the expression level of ST3GAL1 was significantly upregulated in the low-score group. Immunosuppressive cell infiltration represents a potential immune escape mechanism in the low-score group, and targeting these cells may be an effective immunotherapeutic strategy.
Generally, the cytotoxic activity of NK cells is inhibited by the binding of killer-cell immunoglobulin-like receptors to MHC-I molecules on normal cells. In contrast, tumor cells typically lack MHC-I and are thus vulnerable to NK-mediated killing (20, 39, 60). Surprisingly, the high-score group with low MHC-I activity had a low NK cell infiltration and cytotoxicity level, indicating a potent NK cell suppression in this group. The imbalance of NK cell activation and inhibitory signaling is a critical mechanism of immune escape in neuroblastoma (20, 39). We found ligands for NK cell activating receptors, including MICA, MICB, and ULBP1 were downregulated in the high-score group. The inhibition of NK cell activation may be the potential reason for escaping NK cells in the high-score group. Additionally, previous research showed that overexpression of B7-H3 molecules could inhibit NK cell cytotoxicity, and we found that the high-score group exhibited increased B7-H3 expression (60). The escaping NK cell significantly contributes to immunosuppression in the high-score group. Fortunately, the anti-GD2 antibody Dinutuximab could restore the NK cell balance and promote NK cell cytotoxicity against neuroblastoma potently, highlighting the need for anti-GD2 immunotherapy in the high-score group (20).
The presence of immune checkpoints is a critical mechanism by which cancers escape the immune system, and immune checkpoint inhibitors have made breakthroughs in adult malignancies. However, the efficacy of ICIs in neuroblastoma is not satisfactory, and identifying individuals who respond to immunotherapy is essential. We found that PD-L1 was significantly overexpressed in the low-score group, indicating that the low-score group may benefit from ICIs (20). The ganglioside-related risk signature could effectively predict immunotherapy response in the ICI-treated cohort, with better response to immunotherapy in the low-score group. Therefore, ICIs are potential options for immunotherapy in the low-score group, further verification by large-scale and multi-center investigations are required in neuroblastoma.
There are some limitations to this study. There is a shortage of molecular sequencing data for anti-GD2 therapy in neuroblastoma, and the risk signature for predicting anti-GD2 antibody immunotherapy response requires additional confirmation. It is encouraging to note that the anti-GD2 immunotherapy for neuroblastoma has been introduced in China, and a relevant clinical trial is in progress in our institution (CTR20221154). Furthermore, since this study is based on retrospective data, it requires further validation from a comprehensive perspective research. The treatment strategies that we advocated for the different subgroups in the risk signature, including chemotherapy and immunotherapy, should be assessed in prospective clinical trials. Additionally, thorough laboratory investigations are necessary to elucidate the comprehensive biological functions of signature genes.
5. Conclusion
In this work, we developed a novel ganglioside-related risk signature that enabled precise prognostic prediction of neuroblastoma. Additionally, the risk signature identified distinct immune landscapes and immune evasion strategies between risk groups and could be used to predict immunotherapy response. This work emphasizes the critical role of gangliosides in the prognosis and immune microenvironment of neuroblastoma, which may inform clinical evaluation and therapeutic decision-making.
Data availability statement
Publicly available datasets were analyzed in this study. The names of the repositories and accession numbers can be found within the article/Supplementary Materials.
Ethics statement
The studies involving human participants were reviewed and approved by the Ethics Committee of Tianjin Medical University Cancer Institute and Hospital (E20210027). Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.
Author contributions
LH and QZ designed the study. JY, LH, YS, YJ, ZL, BG, JL, YL, and YW analyzed and interpreted data. JY drafted the manuscript, and major revised by LH and QZ. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by grants from the National Key Research and Development Program of China (2018YFC1313000, 2018YFC1313001) and Tianjin Key Medical Discipline (Specialty) Construction Project (TJYXZDXK-009A).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.1061814/full#supplementary-material
Figure S1 | The workflow of the present study.
Figure S2 | Consensus clustering analysis based on differentially expressed ganglioside-related genes in the GSE49710 dataset. (A-H) Consensus score matrix of samples when k = 2-9. (I) CDF of the consensus matrix for each k (indicated by colors). (J) Relative alterations in the area under CDF curves. (K) Tracking plot for each k.
Figure S3 | Consensus clustering analysis based on differentially expressed ganglioside-related genes in the TARGET dataset. (A-H) Consensus score matrix of samples when k = 2-9. (I) CDF of the consensus matrix for each k (indicated by colors). (J) Relative alterations in the area under CDF curves. (K) Tracking plot for each k. (L) The t-SNE analysis revealed a clear distinction between the two clusters in the TARGET cohort. (M) Kaplan–Meier curves of overall survival in TARGET cohort between different ganglioside-related clusters.
Figure S4 | Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis of differentially expressed genes between ganglioside-related clusters A and B (indicated by colors).
Figure S5 | Validation of the ganglioside-related risk signature. (A) The distribution of the risk score (top) and survival time (bottom) in the E-MTAB-8248 cohort. (B) The receiver operating characteristic (ROC) curves of the risk signature for 3-year and 5-year OS prediction in the E-MTAB-8248 cohort. (C) The ROC curves of clinical prognostic factors for 3-year OS prediction in the E-MTAB-8248 cohort. (D) Kaplan–Meier curves of overall survival (OS) in the E-MTAB-8248 cohort between different risk groups. (E) The distribution of the risk score (top) and survival time (bottom) in the TARGET cohort. (F) The ROC curves of the risk signature for 3-year and 5-year OS prediction in the TARGET cohort. (G) The ROC curves of clinical prognostic factors for 3-year OS prediction in the TARGET cohort. (H) Kaplan–Meier curves of overall survival (OS) in the TARGET cohort between different risk groups. (I-K) The principal component analysis (PCA) of the risk signature in GSE49710 (I), E-MTAB-8248 (J), and TARGET (K) cohorts.
Figure S6 | The prognostic value of the ganglioside-related risk signature in samples in COG high-risk group or older than 18 months. (A-B) Kaplan–Meier curves of samples in COG high-risk group in GSE49710 (A) and TARGET (B) cohorts. (C-D) Kaplan–Meier curves of samples older than 18 months in GSE49710 (C) and TARGET (D) cohorts.
Figure S7 | Correlation analysis between clinical characteristics and the risk signature in the Tianjin cohort. (A-E) The comparison of risk scores between samples with different clinical characteristics in the Tianjin cohort, including INRGSS stage (A), MYCN status (B), COG risk groups (C), pathology (D), and INSS stage (E).
Figure S8 | Independent prognosis analysis of the risk signature and external validation of the nomogram. (A-B) The univariate and multivariate Cox regression analysis in E-MTAB-8248 (A) and TARGET (B) cohorts. (C) Calibration curves of the nomogram in the prediction of 3-year and 5-year overall survival (OS) in the TARGET cohort.
Figure S9 | Gene set enrichment analysis (GSEA) of the risk signature. The GSEA in the high-score (A) and low-score (B) groups.
References
1. Kholodenko IV, Kalinovsky DV, Doronin II, Deyev SM, Kholodenko RV. Neuroblastoma origin and therapeutic targets for immunotherapy. J Immunol Res (2018) 2018:7394268. doi: 10.1155/2018/7394268
2. Johnsen JI, Dyberg C, Wickström M. Neuroblastoma-a neural crest derived embryonal malignancy. Front Mol Neurosci (2019) 12:9. doi: 10.3389/fnmol.2019.00009
3. Zafar A, Wang W, Liu G, Wang X, Xian W, McKeon F, et al. Molecular targeting therapies for neuroblastoma: Progress and challenges. Med Res Rev (2021) 41(2):961–1021. doi: 10.1002/med.21750
4. Irwin MS, Naranjo A, Zhang FF, Cohn SL, London WB, Gastier-Foster JM, et al. Revised neuroblastoma risk classification system: A report from the children's oncology group. J Clin Oncol (2021) 39(29):3229–41. doi: 10.1200/jco.21.00278
5. Groux-Degroote S, Guérardel Y, Delannoy P. Gangliosides: Structures, biosynthesis, analysis, and roles in cancer. Chembiochem (2017) 18(13):1146–54. doi: 10.1002/cbic.201600705
6. Groux-Degroote S, Delannoy P. Cancer-associated glycosphingolipids as tumor markers and targets for cancer immunotherapy. Int J Mol Sci (2021) 22(11):6145. doi: 10.3390/ijms22116145
7. Liu J, Zheng X, Pang X, Li L, Wang J, Yang C, et al. Ganglioside Gd3 synthase (Gd3s), a novel cancer drug target. Acta Pharm Sin B (2018) 8(5):713–20. doi: 10.1016/j.apsb.2018.07.009
8. Piazzesi A, Afsar SY, van Echten-Deckert G. Sphingolipid metabolism in the development and progression of cancer: One cancer's help is another's hindrance. Mol Oncol (2021) 15(12):3256–79. doi: 10.1002/1878-0261.13063
9. Wang H, Isaji T, Satoh M, Li D, Arai Y, Gu J. Antitumor effects of exogenous ganglioside Gm3 on bladder cancer in an orthotopic cancer model. Urology (2013) 81(1):210.e11–5. doi: 10.1016/j.urology.2012.08.015
10. Sasaki N, Hirabayashi K, Michishita M, Takahashi K, Hasegawa F, Gomi F, et al. Ganglioside Gm2, highly expressed in the mia paca-2 pancreatic ductal adenocarcinoma cell line, is correlated with growth, invasion, and advanced stage. Sci Rep (2019) 9(1):19369. doi: 10.1038/s41598-019-55867-4
11. Suzuki M, Nagane M, Kato K, Yamauchi A, Shimizu T, Yamashita H, et al. Endothelial ganglioside Gm3 regulates angiogenesis in solid tumors. Biochem Biophys Res Commun (2021) 569:10–6. doi: 10.1016/j.bbrc.2021.06.063
12. Ramos RI, Bustos MA, Wu J, Jones P, Chang SC, Kiyohara E, et al. Upregulation of cell surface Gd3 ganglioside phenotype is associated with human melanoma brain metastasis. Mol Oncol (2020) 14(8):1760–78. doi: 10.1002/1878-0261.12702
13. Cavdarli S, Groux-Degroote S, Delannoy P. Gangliosides: The double-edge sword of neuro-ectodermal derived tumors. Biomolecules (2019) 9(8):311. doi: 10.3390/biom9080311
14. Schengrund CL. Gangliosides and neuroblastomas. Int J Mol Sci (2020) 21(15):5313. doi: 10.3390/ijms21155313
15. Wang J, Yu RK. Interaction of ganglioside Gd3 with an egf receptor sustains the self-renewal ability of mouse neural stem cells in vitro. Proc Natl Acad Sci U.S.A. (2013) 110(47):19137–42. doi: 10.1073/pnas.1307224110
16. Yeh SC, Wang PY, Lou YW, Khoo KH, Hsiao M, Hsu TL, et al. Glycolipid Gd3 and Gd3 synthase are key drivers for glioblastoma stem cells and tumorigenicity. Proc Natl Acad Sci U.S.A. (2016) 113(20):5592–7. doi: 10.1073/pnas.1604721113
17. Shibuya H, Hamamura K, Hotta H, Matsumoto Y, Nishida Y, Hattori H, et al. Enhancement of malignant properties of human osteosarcoma cells with disialyl gangliosides Gd2/Gd3. Cancer Sci (2012) 103(9):1656–64. doi: 10.1111/j.1349-7006.2012.02344.x
18. Sasaki N, Toyoda M, Ishiwata T. Gangliosides as signaling regulators in cancer. Int J Mol Sci (2021) 22(10):5076. doi: 10.3390/ijms22105076
19. Yoshida S, Fukumoto S, Kawaguchi H, Sato S, Ueda R, Furukawa K. Ganglioside G(D2) in small cell lung cancer cell lines: Enhancement of cell proliferation and mediation of apoptosis. Cancer Res (2001) 61(10):4244–52.
20. Wienke J, Dierselhuis MP, Tytgat GAM, Künkele A, Nierkens S, Molenaar JJ. The immune landscape of neuroblastoma: Challenges and opportunities for novel therapeutic strategies in pediatric oncology. Eur J Cancer (2021) 144:123–50. doi: 10.1016/j.ejca.2020.11.014
21. Dillinger B, Ahmadi-Erber S, Lau M, Hoelzl MA, Erhart F, Juergens B, et al. Ifn-Γ and tumor gangliosides: Implications for the tumor microenvironment. Cell Immunol (2018) 325:33–40. doi: 10.1016/j.cellimm.2018.01.014
22. Zhang W, Yu Y, Hertwig F, Thierry-Mieg J, Zhang W, Thierry-Mieg D, et al. Comparison of rna-seq and microarray-based models for clinical endpoint prediction. Genome Biol (2015) 16(1):133. doi: 10.1186/s13059-015-0694-1
23. Roderwieser A, Sand F, Walter E, Fischer J, Gecht J, Bartenhagen C, et al. Telomerase is a prognostic marker of poor outcome and a therapeutic target in neuroblastoma. JCO Precis Oncol (2019) 3:1–20. doi: 10.1200/po.19.00072
24. Kang W, Hu J, Zhao Q, Song F. Identification of an autophagy-related risk signature correlates with immunophenotype and predicts immune checkpoint blockade efficacy of neuroblastoma. Front Cell Dev Biol (2021) 9:731380. doi: 10.3389/fcell.2021.731380
25. Wilkerson MD, Hayes DN. Consensusclusterplus: A class discovery tool with confidence assessments and item tracking. Bioinformatics (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170
26. Yu G, Wang LG, Han Y, He QY. Clusterprofiler: An r package for comparing biological themes among gene clusters. Omics (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118
27. Wang S, He Z, Wang X, Li H, Liu XS. Antigen presentation and tumor immunogenicity in cancer immunotherapy response prediction. Elife (2019) 8:e49020. doi: 10.7554/eLife.49020
28. Liu Z, Guo Y, Yang X, Chen C, Fan D, Wu X, et al. Immune landscape refines the classification of colorectal cancer with heterogeneous prognosis, tumor microenvironment and distinct sensitivity to frontline therapies. Front Cell Dev Biol (2021) 9:784199. doi: 10.3389/fcell.2021.784199
29. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019
30. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612
31. Bhattacharya S, Andorf S, Gomes L, Dunn P, Schaefer H, Pontius J, et al. Immport: Disseminating data to the public for the future of immunology. Immunol Res (2014) 58(2-3):234–9. doi: 10.1007/s12026-014-8516-1
32. Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, et al. Machine learning identifies stemness features associated with oncogenic dedifferentiation. Cell (2018) 173(2):338–54.e15. doi: 10.1016/j.cell.2018.03.034
33. Daily K, Ho Sui SJ, Schriml LM, Dexheimer PJ, Salomonis N, Schroll R, et al. Molecular, phenotypic, and sample-associated data to describe pluripotent stem cell lines and derivatives. Sci Data (2017) 4:170030. doi: 10.1038/sdata.2017.30
34. Salomonis N, Dexheimer PJ, Omberg L, Schroll R, Bush S, Huo J, et al. Integrated genomic analysis of diverse induced pluripotent stem cells from the progenitor cell biology consortium. Stem Cell Rep (2016) 7(1):110–25. doi: 10.1016/j.stemcr.2016.05.006
35. Geeleher P, Cox N, Huang RS. Prrophetic: An r package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS One (2014) 9(9):e107468. doi: 10.1371/journal.pone.0107468
36. Kim D, Langmead B, Salzberg SL. Hisat: A fast spliced aligner with low memory requirements. Nat Methods (2015) 12(4):357–60. doi: 10.1038/nmeth.3317
37. Liao Y, Smyth GK, Shi W. Featurecounts: An efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics (2014) 30(7):923–30. doi: 10.1093/bioinformatics/btt656
38. Cheung IY, Vickers A, Cheung NK. Sialyltransferase stx (St8siaii): A novel molecular marker of metastatic neuroblastoma. Int J Cancer (2006) 119(1):152–6. doi: 10.1002/ijc.21789
39. Vanichapol T, Chutipongtanate S, Anurathapan U, Hongeng S. Immune escape mechanisms and future prospects for immunotherapy in neuroblastoma. BioMed Res Int (2018) 2018:1812535. doi: 10.1155/2018/1812535
40. Reid GS, Shan X, Coughlin CM, Lassoued W, Pawel BR, Wexler LH, et al. Interferon-Gamma-Dependent infiltration of human T cells into neuroblastoma tumors in vivo. Clin Cancer Res (2009) 15(21):6602–8. doi: 10.1158/1078-0432.Ccr-09-0829
41. Chowdhury F, Dunn S, Mitchell S, Mellows T, Ashton-Key M, Gray JC. Pd-L1 and Cd8+Pd1+ lymphocytes exist as targets in the pediatric tumor microenvironment for immunomodulatory therapy. OncoImmunology (2015) 4(10):e1029701. doi: 10.1080/2162402X.2015.1029701
42. Kawasaki Y, Ito A, Withers DA, Taima T, Kakoi N, Saito S, et al. Ganglioside Dsgb5, preferred ligand for siglec-7, inhibits nk cell cytotoxicity against renal cell carcinoma cells. Glycobiology (2010) 20(11):1373–9. doi: 10.1093/glycob/cwq116
43. Nicoll G, Avril T, Lock K, Furukawa K, Bovin N, Crocker PR. Ganglioside Gd3 expression on target cells can modulate nk cell cytotoxicity Via siglec-7-Dependent and -independent mechanisms. Eur J Immunol (2003) 33(6):1642–8. doi: 10.1002/eji.200323693
44. Morandi F, Sabatini F, Podestà M, Airoldi I. Immunotherapeutic strategies for neuroblastoma: Present, past and future. Vaccines (Basel) (2021) 9(1):43. doi: 10.3390/vaccines9010043
45. Morad G, Helmink BA, Sharma P, Wargo JA. Hallmarks of response, resistance, and toxicity to immune checkpoint blockade. Cell (2021) 184(21):5309–37. doi: 10.1016/j.cell.2021.09.020
46. Liang YJ, Wang CY, Wang IA, Chen YW, Li LT, Lin CY, et al. Interaction of glycosphingolipids Gd3 and Gd2 with growth factor receptors maintains breast cancer stem cell phenotype. Oncotarget (2017) 8(29):47454–73. doi: 10.18632/oncotarget.17665
47. Jiang L, Hao Y, Shao C, Wu Q, Prager BC, Gimple RC, et al. Adar1-mediated rna editing links ganglioside catabolism to glioblastoma stem cell maintenance. J Clin Invest (2022) 132(6):e143397. doi: 10.1172/jci143397
48. Battula VL, Shi Y, Evans KW, Wang RY, Spaeth EL, Jacamo RO, et al. Ganglioside Gd2 identifies breast cancer stem cells and promotes tumorigenesis. J Clin Invest (2012) 122(6):2066–78. doi: 10.1172/jci59735
49. Sakai H, Tanaka Y, Tanaka M, Ban N, Yamada K, Matsumura Y, et al. Abca2 deficiency results in abnormal sphingolipid metabolism in mouse brain. J Biol Chem (2007) 282(27):19692–9. doi: 10.1074/jbc.M611056200
50. Aberuyi N, Rahgozar S, Moafi A. The role of atp-binding cassette transporter A2 in childhood acute lymphoblastic leukemia multidrug resistance. Iran J Ped Hematol Oncol (2014) 4(3):118–26.
51. Zhang M, Liu Y, Kong D. Identifying biomolecules and constructing a prognostic risk prediction model for recurrence in osteosarcoma. J Bone Oncol (2021) 26:100331. doi: 10.1016/j.jbo.2020.100331
52. Sha Y, Han L, Sun B, Zhao Q. Identification of a glycosyltransferase signature for predicting prognosis and immune microenvironment in neuroblastoma. Front Cell Dev Biol (2021) 9:769580. doi: 10.3389/fcell.2021.769580
53. Silvestri I, Testa F, Zappasodi R, Cairo CW, Zhang Y, Lupo B, et al. Sialidase Neu4 is involved in glioblastoma stem cell survival. Cell Death Dis (2014) 5(8):e1381. doi: 10.1038/cddis.2014.349
54. Zhang X, Yang X, Chen M, Zheng S, Li J, Lin S, et al. St3gal3 confers Paclitaxel−Mediated chemoresistance in ovarian cancer cells by attenuating Caspase−8/3 signaling. Mol Med Rep (2019) 20(5):4499–506. doi: 10.3892/mmr.2019.10712
55. Hugonnet M, Singh P, Haas Q, von Gunten S. The distinct roles of sialyltransferases in cancer biology and onco-immunology. Front Immunol (2021) 12:799861. doi: 10.3389/fimmu.2021.799861
56. Wu X, Zhao J, Ruan Y, Sun L, Xu C, Jiang H. Sialyltransferase St3gal1 promotes cell migration, invasion, and tgf-β1-Induced emt and confers paclitaxel resistance in ovarian cancer. Cell Death Dis (2018) 9(11):1102. doi: 10.1038/s41419-018-1101-0
57. Guerrero PE, Miró L, Wong BS, Massaguer A, Martínez-Bosch N, Llorens R, et al. Knockdown of α2,3-sialyltransferases impairs pancreatic cancer cell migration, invasion and e-Selectin-Dependent adhesion. Int J Mol Sci (2020) 21(17):6239. doi: 10.3390/ijms21176239
58. Spel L, Nieuwenhuis J, Haarsma R, Stickel E, Bleijerveld OB, Altelaar M, et al. Nedd4-binding protein 1 and Tnfaip3-interacting protein 1 control mhc-1 display in neuroblastoma. Cancer Res (2018) 78(23):6621–31. doi: 10.1158/0008-5472.Can-18-0545
59. Durgeau A, Virk Y, Corgnac S, Mami-Chouaib F. Recent advances in targeting Cd8 T-cell immunity for more effective cancer immunotherapy. Front Immunol (2018) 9:14. doi: 10.3389/fimmu.2018.00014
Keywords: neuroblastoma, gangliosides, prognosis, immune landscape, dedifferentiation
Citation: Yang J, Han L, Sha Y, Jin Y, Li Z, Gong B, Li J, Liu Y, Wang Y and Zhao Q (2022) A novel ganglioside-related risk signature can reveal the distinct immune landscape of neuroblastoma and predict the immunotherapeutic response. Front. Immunol. 13:1061814. doi: 10.3389/fimmu.2022.1061814
Received: 05 October 2022; Accepted: 05 December 2022;
Published: 20 December 2022.
Edited by:
Ping Zheng, The University of Melbourne, AustraliaReviewed by:
Valeria I Segatori, National University of Quilmes, ArgentinaShiyu Mao, Tongji University, China
Copyright © 2022 Yang, Han, Sha, Jin, Li, Gong, Li, Liu, Wang and Zhao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Qiang Zhao, emhhb3FpYW5nQHRqbXVjaC5jb20=
†These authors contributed equally to this work and share first authorship