Identification of a Glycosyltransferase Signature for Predicting Prognosis and Immune Microenvironment in Neuroblastoma

Neuroblastoma (NB) is one of the most common solid tumors in children. Glycosyltransferases (GTs) play a crucial role in tumor development and immune escape and have been used as prognostic biomarkers in various tumors. However, the biological functions and prognostic significance of GTs in NB remain poorly understood. The expression data from Gene Expression Omnibus (GEO) and Therapeutically Applicable Research to Generate Effective Treatments (TARGET) were collected as training and testing data. Based on a progression status, differentially expressed GTs were identified. We constructed a GTscore through support vector machine, least absolute shrinkage and selection operator, and Cox regression in NB, which included four prognostic GTs and was an independent prognostic risk factor for NB. Patients in the high GTscore group had an older age, MYCN amplification, advanced International Neuroblastoma Staging System stage, and high risk. Samples with high GTscores revealed high disialoganglioside (GD2) and neuron-specific enolase expression levels. In addition, a lack of immune cell infiltration was observed in the high GTscore group. This GTscore was also associated with the expression of chemokines (CCL2, CXCL9, and CXCL10) and immune checkpoint genes (cytotoxic T-lymphocyte–associated protein 4, granzyme H, and granzyme K). A low GTscore was also linked to an enhanced response to anti–PD-1 immunotherapy in melanoma patients, and one type of tumor was also derived from neuroectodermal cells such as NB. In conclusion, the constructed GTscore revealed the relationship between GT expression and the NB outcome, GD2 phenotype, and immune infiltration and provided novel clues for the prediction of prognosis and immunotherapy response in NB.

Introduction

Neuroblastoma (NB) is a malignant childhood solid tumor of neuroectodermal cells derived from the sympathetic nervous system and the primary cause of pediatric cancer-related mortality (Maris et al., 2007). Despite advances in multimodality therapy, the prognosis of high-risk NB patients remains poor. Age, tumor stage, mitosis-karyorrhexis index (MKI), MYCN amplification, and ALK expression are used to predict the clinical outcomes of NB (Perel et al., 2004; Chang et al., 2020; Sokol et al., 2020). However, there are limitations in the accuracy and specificity of each index. Therefore, it is urgent to explore new prognostic markers.

Glycosyltransferases (GTs) are a series of enzymes that catalyze the formation of glycosidic linkages to diverse acceptor molecules, such as proteins, lipids, hormones, and oligosaccharides. The altered expression of GTs leading to aberrant glycosylation patterns has been found in various cancers and is considered as a biomarker of cancer, such as small cell lung cancer (Yoshida et al., 2001), NB (Berois et al., 2013), melanoma (Hu et al., 2019), and triple-negative breast cancer (TNBC) (Narvaez et al., 2020). Moreover, the glycosylation products catalyzed by various GTs also play an important role in the biological function of cancer cells, including proliferation, migration, metastasis, immune escape, and apoptosis (Cavdarli et al., 2019). A large number of studies have been performed to investigate the functions of GTs in various types of tumors (Huang et al., 2015; Schultz et al., 2016; Dong et al., 2017; Bhat et al., 2018). Unfortunately, to date, there is still a lack of an in-depth analysis of the role of these GTs in NB. With large amounts of molecular sequence data accumulating gradually, bioinformatics based on the massive microarray datasets provides a promising strategy for these studies.

Disialoganglioside (GD2), which is synthesized by GTs, is highly expressed on NB tumor cells, with absent or weak expression in normal tissues (Balis et al., 2020). In addition to NB, GD2 is also overexpressed in other neuroectodermal tumors and sarcomas, such as melanoma (Hamilton et al., 1993), brain tumors (Mennel et al., 2000), osteosarcoma (Poon et al., 2015), and Ewing sarcoma (Fisher et al., 2016), as well as several types of breast cancer (Orsi et al., 2017) and small cell lung cancer (Cheresh et al., 1986). GD2 has many biological functions, involving tumor cell proliferation, migration, invasion, and apoptosis (Liu et al., 2014; Mansoori et al., 2019). Schulz et al. showed that soluble GD2 (sGD2) levels in NB patients were associated with disease progression (Schulz et al., 1984). Furthermore, tumor-derived GD2 inhibits T-cell proliferation (Ladisch et al., 1992), and therefore may promote tumor immune evasion. Given its high prevalence in NB tissues compared with normal tissues, GD2 is used as a target for cancer immunotherapy. As GD2-positive tumors are more responsive to such targeted therapy, GD2 has been recommended as a therapeutic monitoring tool in NB patients treated with GD2-targeted immunotherapy (Terzic et al., 2018). Neuron-specific enolase (NSE) is currently considered a useful marker in the diagnosis, prognosis, and follow-up of related neuroendocrine tumors (Zeltzer et al., 1983; Shine et al., 1990; Giovanella et al., 1995). In view of the biological importance of GD2 and NSE, it is imperative to detect the expression levels of GD2 and NSE in cancer cells.

In this study, the expression profiles of GT genes in the Gene Expression Omnibus (GEO) database (GSE49710 dataset) were downloaded and analyzed to identify differentially expressed GTs in NB. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses revealed that genes were mainly involved in the glycosphingolipid biosynthetic process. We systemically assessed the identified differentially expressed GTs, especially the genes involved in the GD2 synthesis pathway. Then, using the least absolute shrinkage and selection operator (LASSO), Cox, and support vector machine (SVM) algorithm, we developed a GTscore for predicting the prognosis of NB patients via the GSE49710 dataset and verified the feasibility of the GTscore through the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) dataset. The GTscore has a relation with clinicopathological risk factors and could be used as an independent risk factor. Furthermore, our GTscore correlated with immune cell infiltration, GD2 molecular phenotypes, and NSE expression levels. Importantly, the GTscore could accurately predict the response to immunotherapy in melanoma, another type of tumor derived from neuroectodermal cells. In summary, we built a 4-GT signature model that could be used as a predictor of NB survival, responsiveness to immunotherapy in melanoma, and GD2 molecular phenotypes of tumor cells.

Materials and Methods

Data Sources and Data Preprocessing

The gene expression data and clinical information of GSE49710 (Zhang et al., 2015) were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). A total of 498 NB samples including RNA-seq and clinical information (sex, age, International Neuroblastoma Staging System [INSS] stage, MYCN amplification status, and progression) were used as the training set. The GSE112447, GSE90689, and GSE142293 datasets were also downloaded for GD2 molecular phenotype analysis. Probes were converted to corresponding gene symbols according to the “clusterProfiler” R package (Yu et al., 2012) with the platforms GPL16876 and GPL6480 to obtain a gene expression matrix. The gene expression profile and clinical data (sex, age, INSS stage, MYCN amplification status, MKI, grade, COG risk, relapse, and progression) of TARGET-NB were downloaded from the TARGET database (154 patients) (https://ocg.cancer.gov/programs/target), which were used for validation of differentially expressed genes. Gene sequencing data and immunotherapy information of The Cancer Genome Atlas Skin Cutaneous Melanoma (TCGA-SKCM) were downloaded from TCGA database (https://portal.gdc.cancer.gov/). Ensembl IDs were converted to corresponding gene symbols with the clusterProfiler R package. Then, these gene expression matrices were quantile standardized and log2 transformed for subsequent analysis.

Screening Differentially Expressed GTs

The processed gene expression matrix of GSE49710 was subjected to differentially expressed GT gene (DEG) analysis to screen genes that were differentially expressed between the progression and nonprogression groups. The “limma” (Ritchie et al., 2015) R package was utilized, and genes were identified as significant DEGs if |log2FC| value >0.5 and false discovery rate (FDR) <0.05. Then, univariable Cox regression analysis was applied to identify DEGs related to the prognosis of NB patients.

Functional Enrichment Analysis

In this study, the clusterProfiler package was used to complete the enrichment analysis with DEGs, aiming to provide a primary understanding of the biological function of these genes. p < 0.05 and FDR < 25% were identified as significant.

Protein-Protein Interaction Network Construction

The STRING (version 11.0, https://string-db.org) database was used to construct protein–protein interaction (PPI) network (Szklarczyk et al., 2015). Then, the network was visualized via Cytoscape software (version 3.5.1).

Feature Selection and Identification of GTscore

The LASSO is a shrinkage estimation method whose basic premise is to minimize the residual sum of squares under the constraint by restricting coefficients of some features to zero and further obtain shrinkage subsets (Tibshirani, 1997). In this study, a LASSO Cox regression model was applied to identify DEGs associated with progression-free survival (PFS). Tenfold cross-validation for tuning parameter selection was performed, and the partial likelihood deviance met the minimum criteria with the “glmnet” R package.

The SVM method is used to find the best variable combination by subtracting the inputted feature vector (Huang et al., 2018). In the current study, the gene expression matrix of 498 pediatric NB samples was subjected to SVM operation and randomly divided into a training group and verification group with the “glmnet” and “e1071” R packages. The ganglioside synthesis–related DEGs were combined in various combinations to establish the optimal predictive classification model, and the smallest classification error and highest accuracy of the prediction model of the training group were externally verified to obtain the optimal combination.

The advantage of these algorithms is that they can retain subset shrinkage, so they are more suitable for processing multicollinearity data. We then obtained the coefficients of screened DEGs to construct a GTscore. The risk score for the individual samples was calculated as follows:

$\text{GTscore }=\text{ β}1\text{∗X}1+\text{β}2\text{∗X}2+\dots \text{βi∗Xi }\dots +\text{ βn∗Xn}$

where β_i represents the correlation coefficient for gene i, and X_i represents the expression of gene i.

According to the median cutoff value of the GTscore, the patients were divided into high and low GTscore groups. The log-rank test was applied to compare the overall survival (OS), event-free survival (EFS), and recurrence-free survival (RFS) differences between the two groups. In addition, univariate Cox regression and multivariable Cox regression were used to assess the independent prognostic value of the GTscore.

Immune Cell Infiltration Evaluation

It is generally believed that the degree of immune cell infiltration is associated with tumor occurrence, development, treatment, and clinical prognosis (Whiteside, 2008; Asgharzadeh et al., 2012; Wienke et al., 2021). Studies have found that GT genes are related to tumor immune escape (Wang et al., 2019; Liu et al., 2020). Accordingly, we explored the relationship between the GTscore and immune cell infiltration. With previously published marker gene sets representing diverse immune cell types (Bindea et al., 2013), the relative abundance of immune cell types for a single sample was quantified using the R package “GSVA” (Hanzelmann et al., 2013) and single-sample gene set enrichment analysis (ssGSEA) function under standard settings. To further explore the inflammatory status, we estimated the expression levels of a series of chemokines in the high and low GTscore groups. In addition, the relationships between the expression levels of immune checkpoint genes such as cytotoxic T-lymphocyte–associated protein 4 (CTLA-4) and granzymes H and K (GZMH, GZMK) and the GTscore were evaluated.

Correlation of the GTscore With GD2 Molecular Phenotypes and NSE Expression Levels

GD2 and NSE were the significant markers for neuroendocrine tumors, including NB, melanoma, and small cell lung cancer. Because of the importance of GD2 and NSE, we assessed the correlation between GTscore and GD2 molecular phenotypes and NSE expression levels.

To evaluate the relationship between GTscore and GD2 molecular phenotypes, we compared the GTscore between putatively GD2 (+) (GSE112447) and GD2 (−) NB samples (GSE90689), GD2 (++) cell lines, and GD2 (+/−) cell lines (GSE142293) (Sorokin et al., 2020), GD2 (−) acute lymphoblastic leukemia (ALL) (Macher et al., 1982; Pick et al., 1986) and GD2 (+) NB (Wu et al., 1986; Terzic et al., 2018) (TARGET database), GD2 (−) non-TNBC, and GD2 (+) TNBC samples (Orsi et al., 2017) (TCGA database).

Furthermore, Spearman correlation analysis was used to analyze the relationship between GTscore and NSE expression levels in NB (Zhang et al., 2015) (GSE49710), melanoma (Bakhoum and Esmaeli, 2019) (TCGA database), and small cell lung cancer (TCGA database).

Estimation of Immunotherapeutic Response

Tumor immunotherapy using immune checkpoint inhibitors (ICIs) is a promising treatment for solid tumors. GTs play an important role in cancer cell immune escape (Galvan et al., 2000; Wang et al., 2019). Thus, we explored the correlation between GTscore and tumor response to ICIs in SKCM, a type of neuroendocrine tumor. Receiver operating characteristic (ROC) curves were applied to investigate the predictive ability of the GTscore for immunotherapeutic response. The Kaplan–Meier method was used to analyze 1- and 3-year PFS stratified by GTscore.

Statistical Analysis

Student t test, analysis of variance, Mann–Whitney Wilcoxon test, and Kruskal–Wallis test were used to compare continuous variables and ordered categorical variables. Pearson χ² test and Fisher exact test were used to compare unordered categorical variables. Survival curves were constructed using the Kaplan–Meier method and were compared using the log-rank test and “survminer” R package. The optimal cutoff point for continuous variables was determined by the “surv_cutpoint” function in the survminer R package. A Cox proportional hazard regression model was performed to calculate hazard ratios (HRs) and 95% confidence intervals (CIs). Multiple testing corrections were performed using the Benjamini–Hochberg procedure and produced FDRs. All the tests were two-sided, and p < 0.05 or FDR < 25% was regarded as indicating statistical significance unless otherwise stated. All analyses were performed using R package version 3.4.2 (https://cran.r-project.org/).

Results

Identification of Differentially GT Expressed Genes Associated With NB Progression and Patient Survival

The workflow is shown in Figure 1. First, 213 GT genes were collected from MSigDB (Molecular Signatures Database) (https://www.gsea-msigdb.org/gsea/msigdb). We then explored which genes were differentially expressed between progression and nonprogression NB samples based on the dataset GSE49710, the cohort of 498 NB samples with gene expression data, and clinicopathological information (Table 1). As shown in Figures 2A–C, we found that 64 GT genes were significantly related to NB progression and survival. GO analysis results showed that DEGs were mostly enriched in glycosphingolipid biosynthetic process, protein glycosylation, and galactosyltransferase activity (Figure 2D). KEGG analysis results suggested that these DEGs were mainly involved in glycosphingolipid biosynthesis and glycan biosynthesis (Figure 2E). The results above suggested that the glycosphingolipid biosynthesis biological process played an important role in the progression of NB. Anti-GD2 immunotherapy showed favorable clinical efficacy for high-risk NB (Rovigatti, 2020). GD2 is synthesized under the catalysis of GTs (Groux-Degroote et al., 2017), which are part of the DEGs. Thus, we focused on the 10 GD2 synthesis–related GT genes in subsequent analysis.

FIGURE 1

FIGURE 1. Flowchart of this study.

TABLE 1

TABLE 1. Clinical pathological parameters of NB patients in the training set.

FIGURE 2

FIGURE 2. DEGs in GSE49710 dataset. (A) Volcano plot shows DEGs between progression and nonprogression NB samples. Red points represent upregulated DEGs. Green points denote downregulated DEGs. Gray points denote genes with no significant differences. (B) Heat map showing the association between 64 DEGs and clinical characteristics annotated on the top. (C) Forest plot of univariate Cox regression results of DEGs. GO (D) and KEGG (E) analyses of DEGs.

Construction and Validation of GT Gene Prognostic Model

Difference analysis and survival analysis illustrated that these GD2 synthesis–related GT genes were significantly related to EFS (Supplementary Figure S1) and OS (Supplementary Figure S2). B4GALT5 and ST3GAL2 were considered as risk factors, whereas the remaining eight genes acted as protective factors. The PPI network and expression correlation analysis showed that they were highly intercorrelated (Supplementary Figure S3). With the SVM algorithm, we screened out the optimal gene combination: ST3GAL2, ST3GAL3, B4GALT5, and B3GALT4 (Figures 3A,B). Next, we conducted LASSO Cox regression analysis and obtained coefficients for each gene (Figures 3C–E). The final GTscore formula was as follows: GTscore = (0.4423) * B4GALT5 expression value + (−0.3630) * B3GALT4 expression value + (0.0636) * ST3GAL2 expression value + (−0.4962) * ST3GAL3 expression value. The ROC curve analyses revealed that the area under the curve (AUC) values at 1-, 3-, and 5-year EFS and OS for the GTscore were higher than those for individual genes (Supplementary Figure S4), indicating that the prognostic role of the GTscore was more accurate than that of a single gene alone.

FIGURE 3

FIGURE 3. Construction of GTscore with machine learning in the training set. The accuracy (A) and the error (B) of the estimate generation for the SVM-RFE algorithm. (C) Ten-time cross-validation was performed for tuning parameter selection through the LASSO regression model. (D) LASSO coefficient profiles of prognostic GTs. (E) Forest plot showing the HR with 95% CI of prognostic GTs.

In the training cohort, Kaplan–Meier survival analysis results showed that the EFS and OS of patients in the high GTscore group were significantly lower than those in the low GTscore group (HR = 2.64, p < 0.001 for EFS, and HR = 6.15, p < 0.001 for OS) (Figures 4A, B). The AUCs of the GTscore to distinguish 1-, 3-, and 5-year EFS were 0.619, 0.677, and 0.661, respectively (Figure 4C), and the AUCs to distinguish 1-, 3-, and 5-year OS were 0.812, 0.785, and 0.779, respectively (Figure 4D), demonstrating great efficacy in distinguishing the prognostic status of NB. A summary of the clinical characteristics of the validation cohort is shown in Table 2. In the validation cohort, Kaplan–Meier survival analysis results showed that the RFS and OS of patients in the high GTscore group were significantly lower than those in the low GTscore group (HR = 1.91, p = 0.026 for RFS, and HR = 1.65, p = 0.041 for OS) (Figures 5A, B). The AUCs to distinguish 1-, 3-, and 5-year RFS were 0.870, 0.732, and 0.724 (Figure 5C), and 0.643, 0.659, and 0.709 for 1-, 3-, and 5-year OS (Figure 5D), suggesting that the GTscore was a robust prognostic model and had good predictive performance in different circumstances.

FIGURE 4

FIGURE 4. Validation of the prognostic value of the GTscore for NB patients in the training set. Kaplan–Meier analysis for EFS (A) and OS (B) of NB patients based on the risk stratification. ROC analysis for 1-, 3-, and 5-year EFS (C) and OS (D) prediction in NB patients. The AUCs at 1, 3, and 5 years were calculated to assess the predictive accuracy.

TABLE 2

TABLE 2. Clinical pathological parameters of NB patients in the validation set.

FIGURE 5

FIGURE 5. Validation of the prognostic value of GTscore for NB patients in the validation set. Kaplan–Meier analysis for RFS (A) and OS (B) of NB patients. ROC analysis for 1-, 3-, and 5-year RFS (C) and OS (D) prediction in NB patients. The AUCs at 1, 3, and 5 years were calculated to assess the predictive accuracy.

The GTscore Is an Independent Prognostic Factor

Currently, clinicopathological factors are still the most important guidelines and are used to predict the prognosis of NB. Therefore, we explored the correlation between GTscore and clinical characteristics. Correlation analysis showed that GTscore was significantly related to age, INSS staging, MYCN amplification, high risk, and disease progression in the training cohort (Figure 6A, Supplementary Figure S5A), and MYCN amplification in the validation cohort (Figure 6B, Supplementary Figure S5B) (all p < 0.05).

FIGURE 6

FIGURE 6. GTscore is related to the clinical characteristics of NB and is an independent prognostic factor. The correlation heat map showing the association of GTscore with clinical characteristics in the training set (A) and validation set (B). Univariate analysis and multivariate analysis of the GTscore, age, sex, INSS stage, MYCN amplification, clinical risk, disease progression, MKI, and OS in the training set (C) and validation set (D).

In the training cohort, univariate Cox regression results showed that patients with older age (>18 months), higher INSS stage, MYCN amplification, higher clinical risk, disease progression, and higher GTscore had a worse prognosis (all p < 0.05) (Figure 6C). Subsequently, the multivariate Cox regression analysis showed that high INSS staging, MYCN amplification, high clinical risk, disease progression, and high GTscore (p < 0.001) were independent risk factors for NB patients (Figure 6C).

In the validation cohort, the correlation analysis showed that patients with high GTscore were older (>18 months), had higher INSS stage and MYCN amplification, had higher MKI (>4%), and had higher clinical risk and disease progression (all p < 0.05) (Figure 6D). Subsequently, these significant indicators were further included in the multivariate Cox regression analysis, and the results showed that disease progression (HR = 3.90, p < 0.001) and high GTscore (HR = 1.35, p = 0.021) were independent risk factors for NB patients (Figure 6D).

Correlation of the GTscore With Immune Cell Infiltration

GTs are closely associated with tumor immune escape (Wang et al., 2019). The tumor microenvironment, which constitutes various immune cell subsets, affects the antitumor effects of immunotherapy (Liu and Joshi, 2020). Therefore, a difference in immune cell infiltration between the two groups was observed. To compare immune cell infiltration between separate GTscore groups, we utilized the ssGSEA approach to explore the relative abundance of immune cell types based on expression profiling data and curated immune gene signatures. Unsupervised hierarchical clustering was applied to categorize the cohort into two infiltration subgroups: a high infiltration group and a low infiltration group. We then investigated the association of immune cell infiltration levels and GTscore. The results showed that most of the 498 NB samples had a low degree of immune cell infiltration (Figure 7A), and the high GTscore group was related to low immune cell infiltration patients (χ² = 38.795, p < 0.001) (Figure 7B). The immune cell network showed that the content proportion of T-cell subgroups was significantly related to the OS of NB patients (p < 0.05) (Figure 7C). We found that the GTscore was negatively correlated with the expression levels of CD56^dim natural killer (NK) cells, CD8⁺ T cells, cytotoxic T cells, dendritic cells (DCs), neutrophils, and T cells but also negatively correlated with the levels of macrophages and regulatory T (Treg) cells (all p < 0.001) (Figure 7D). Furthermore, we found that the expression level of CTLA-4 was significantly positively correlated with the GTscore, whereas that of GZMH and GZMK was significantly negatively correlated with the GTscore (all p < 0.05) (Figure 7E). Considering the role of chemokines and immune checkpoints in immune cell infiltration and activation, we also investigated the association of their expression levels with the GTscore. As shown in Figure 7F, the chemokines CCL2, CXCL9, and CXCL10 were significantly different between the low and high GTscore groups.

FIGURE 7

FIGURE 7. Correlation of GTscore with immune cell infiltration. (A) Heat map of the association between the risk score and the abundance of 24 immune cell types in NB. (B) Sankey plot showing the association between GTscore and immune infiltration. (C) Cellular interaction of the immune cell types in NB. The thickness of the line represents the strength of the correlation. Correlation between GTscore and immune cell infiltration (D), immune checkpoint genes (E), and the expression levels of chemokines (F). ^*p < 0.05, ^***p < 0.001.

The Predictive Value of GTscore for GD2 Phenotypes and NSE Expression Levels in NB and Other Neuroendocrine Tumors

As GD2-positive tumors are more responsive to anti-GD2 targeted immunotherapy, it is imperative to detect the expression levels of GD2 in tumor cells. Therefore, we explored the relationship between GTscore and GD2 phenotypes. In the GSE142293 dataset, the authors analyzed GD2 status of distinct cancer cell lines, including GD2 (−) cell line HOS; GD2 (+) cell lines U2OS, U-373, SH-SY5Y; and GD2 (++) cell lines T98G, IMR-32. We calculated the GTscore for each cell line, and the results showed that the GTscore of the GD2 (++) cell lines was significantly higher than that of the GD2 (+) and GD2 (−) cell lines (p < 0.05) (Figure 8A). We analyzed the difference in GTscore between GD2 (−) ALL and GD2 (+) NB (TARGET database) and GD2 (−) non-TNBC and GD2 (+) TNBC samples (TCGA database). The results showed that GTscores were significantly higher in GD2 (+) samples, corresponding to the GD2 level (p < 0.05) (Figures 8B, D). Moreover, we analyzed the difference in GTscore between flow cytometry–identified 11 GD2 (+) and 24 GD2 (−) NB samples in datasets GSE90689 and GSE112447. The results also demonstrated that the GTscore was significantly higher in GD2 (+) samples (p < 0.05) (Figure 8C).

FIGURE 8

FIGURE 8. Correlation of GTscore with GD2 phenotype and NSE expression levels. (A) Differences in GTscore between the GD2 (−) HOS cell line, GD2 (+) U2OS cell line, U-373 cell line, SH-SY5Y cell line, and GD2 (++) T98G and IMR-32 cell lines. (B) Differences in GTscore between GD2 (−) non-TNBC and GD2 (+) TNBC samples. (C) Differences in GTscore between GD2 (+) and GD2 (−) NB samples. (D) Differences in GTscore between GD2 (−) ALL and GD2 (+) NB samples. Scatterplots showing the correlations between GTscore and NSE expression levels in NB (E), melanoma (F), and small cell lung cancer (G) samples. ^*p < 0.05, ^***p < 0.001, ^****p < 0.0001.

NSE is currently considered a useful marker in the diagnosis and prognosis of neuroendocrine tumors (Braga et al., 2013). Therefore, we investigated the relationship between GTscore and NSE expression levels in NB (GSE49710), melanoma (TCGA database), and small cell lung cancer (TCGA database). The results showed that the GTscore was positively correlated with the expression levels of NSE (all p < 0.001) (Figures 8E–G), suggesting that the GTscore could be used as a marker of neuroendocrine tumors.

GTscore as a Predictive Indicator of Response to Immunotherapy of Neuroendocrine Tumors

ICI therapy has been used for cancer patients (Bagchi et al., 2021). Unfortunately, only a minority of patients can benefit from immunotherapy. As mentioned previously, the GTscore was negatively related to the degree of cytotoxic T-cell infiltration, which was reported to be related to the immunotherapeutic effects of NB patients (Figure 7C). Furthermore, we found that the GTscore was associated with GD2 and NSE, significant markers for neuroendocrine tumors. Because the ICI treatment dataset was not available in NB samples, we chose skin melanoma to investigate whether the GTscore could predict immunotherapeutic benefits in the ICI cohort of skin melanoma (TCGA-SKCM).

According to the evaluation criteria for the efficacy of solid tumors, the results showed that GTscore was negatively related to the clinical efficacy (Figure 9A), and the GTscore of CR samples was significantly lower than that of PD and SD samples (p < 0.05) (Figure 9B). The 1- and 3-year PFS of high GTscore patients was lower than that of low GTscore patients (HR = 3.20 [95% CI = 1.04–9.87, p = 0.038] and HR = 2.40 [CI = 1.05–5.45, p = 0.039], respectively) (Figures 9C, D). The AUC value of the GTscore predicting the benefits of anti–PD-1 immunotherapy was 0.733 (Figure 9E). These results suggested the authentic potency of the GTscore as a biomarker for immunotherapeutic benefit predictions and helped to develop novel treatment strategies for neuroendocrine tumors.

FIGURE 9

FIGURE 9. GTscore as a predictive indicator of the response to immunotherapy of melanoma. (A) Waterfall plot showing the GTscore and corresponding immunotherapy efficacy of individual patients. (B) Boxplot showing the distribution of GTscore in groups with different anti–PD-1 clinical response statuses. GTscore predicting the 1-year (C) and 3-year (D) OS rate of patients receiving anti–PD-1 immunotherapy via Kaplan–Meier curves. (E) The predictive value of GTscore in patients receiving anti–PD-1 immunotherapy via ROC curve analysis. ^*p < 0.05, ^**p < 0.01.

Discussion

NB, the most common childhood extracranial solid tumor, has a poor prognosis, despite comprehensive treatment. Hence, it is necessary to explore the novel prognostic indicators to improve the prognosis of NB patients. Bioinformatics analysis methods provide a new idea for the evaluation of the prognosis of NB patients. A number of models to predict the prognosis of NB, including MYCN-related signatures (Wang et al., 2020), hypoxia gene signatures (Fardin et al., 2010), and immune-based models (Song et al., 2020), have been reported. A large number of studies have shown that GTs are related to cancer progression and can potentially serve as novel biomarkers (Potapenko et al., 2010; Burchell et al., 2018; Tang et al., 2019). GT-based signatures have been constructed to predict the survival in head and neck squamous cell carcinoma (Jin and Yang, 2019) and pancreatic cancer (Gupta et al., 2020). However, the association between GT-related signatures and the prognosis of patients with NB is still unknown. To our knowledge, this study is the first time to identify GTs associated with the prognosis of NB by bioinformatics methods.

In this study, we identified 64 GT DEGs associated with the progression and prognosis of NB patients in the GEO dataset (GSE49710). Moreover, these genes were enriched in the glycosphingolipid biosynthetic process. GD2 is an important glycosphingolipid that has been used as a target for NB immunotherapy (Suzuki and Cheung, 2015). Among them, 10 genes were associated with the GD2 synthesis pathway. Subsequently, based on 10 genes, we identified a GTscore including four GTs using LASSO, Cox regression, and SVM analysis. A novel 4-GT gene signature was then built. Each patient was assigned a score based on GTscore and then assigned to high and low GTscore groups by the optimal cutoff risk scores. In the training set, our GTscore showed a strong ability to predict the prognosis of NB. In addition, the prognostic power of the GTscore was verified in the validation set. Currently, clinicopathological factors are still the most important guidelines and are used to predict the prognosis of NB. Our studies revealed that our GTscore was closely associated with many malignant clinical features. Of note, we found that the GTscore was also an independent prognostic factor. Therefore, our GTscore could serve as a reliable prognostic indicator.

Several studies have shown that GD2 is an important cell-surface molecule involved in cancer biological processes, including tumor cell proliferation, metastasis, invasion, and signal transduction (Mujoo et al., 1987; Iwasawa et al., 2018; Ohmi et al., 2018). GD2 is synthesized from simpler ganglioside molecules by several GTs (Groux-Degroote et al., 2017). Our GTscore comprises four genes, including B4GALT5, B3GALT4, ST3GAL2, and ST3GAL3. Regrettably, few of them were identified as prognostic biomarkers in NB. However, there are still some studies involved in the function of these genes in other tumors. B4GALT5 is a member of the β1,4-galactosyltransferase family. It has been reported that B4GALT5 can suppress apoptosis and enhance cancer cell proliferation (Xu et al., 2002; Jiang and Gu, 2010). B4GALT5 can also mediate the stemness of breast cancer and is considered to be a prognostic biomarker for breast cancer (Tang et al., 2020). Chatterjee et al. reported that B4GALT5 might be a candidate diagnostic and therapeutic biomarker for colorectal cancer (Chatterjee et al., 2019). The B3GALT4 gene belongs to the β-1,3-galactosyltransferase (β3GalT) gene family. It was overexpressed in multiple malignant tumor tissues. Seko et al. had confirmed that B3GALT4 was a promising biomarker for the diagnosis of gynecological cancers (Seko et al., 2009). Zhang et al. also found that B3GALT4 was identified as a prognostic and therapeutic biomarker in colorectal cancer (Zhang et al., 2018). α2,3-Sialyltransferases ST3GAL2 and ST3GAL3 are responsible for sialyl-Lewis(x) synthesis, which is associated with cancer aggressiveness. It has been found that ST3GAL2 is a poor prognostic factor in breast cancer patients treated with chemotherapy (Aloia et al., 2015). Saito et al. also reported that ST3GAL2 activation enhanced the malignant progression of renal cell carcinoma (Saito et al., 2003). Mehta et al. identified that ST3GAL2 and ST3GAL3 overexpression was significantly associated with tumor stage, invasiveness, and metastasis of oral squamous cell carcinoma (Mehta et al., 2020). In pancreatic adenocarcinoma, overexpression of ST3GAL3 correlates with tumor progression (Perez-Garay et al., 2010). These studies supported our GTscore as a prognostic factor for NB patients.

Tumor microenvironment, including stromal cells and immune cells, is closely related to tumor progression in NB (Joshi, 2020). Tumor-infiltrating immune cells are related to disease progression, therapeutic responses, and prognosis (Luen et al., 2017; Ding et al., 2018; Lin et al., 2020). Cancer cells escape from immune surveillance through the imbalance of infiltrating immune cells in the tumor microenvironment (Vanichapol et al., 2018). Consequently, we investigated the variations in immune cell infiltration between the high and low GTscore groups. The results indicated that the numbers of CD8⁺ T cells, cytotoxic T cells, DCs, CD56^dim NK cells, and T cells were significantly lower in the high GTscore group, indicating that our signature may interfere with immune cell infiltration in NB. Gangliosides regulate T-cell activation by affecting the formation of the cell membrane signal complex (Sorice et al., 2004). It has been reported that inhibition of B4Galt5 and B4Galt6 results in higher invariant (i) NK T (iNKT) cell–activating capacities of adipocytes (Rakhshandehroo et al., 2019). It has been reported that overexpression of CD8⁺ T cells or CD56^dim NK cells is associated with favorable prognosis, whereas Treg cells are involved in an unfavorable prognosis (Johnson et al., 2007; Mina et al., 2015; Troschke-Meurer et al., 2019). However, contrary to previous research results, we found that the GTscore was negatively correlated with the levels of macrophages and Treg cells. Thus, further detailed analyses are required to assess the relationships between GTscore and immune cell infiltration. Chemokines play an important role in immune cell infiltration and immunotherapy. Interestingly, CCL2, CXCL9, and CXCL10 were significantly different between the low and high GTscore groups. It has been reported that CCL2 could mediate the migration of iNKT cells to NB in vitro (Metelitsa et al., 2004). It is generally believed that CXCL9 and CXCL10 play an important role in the recruitment and activation of T cells. Importantly, Bocca et al. found that CXCL10 could attract CXCR3⁺ GD2-CAR T cells in vitro (Bocca et al., 2017). In addition, several immune checkpoints including GZMH, GZMK, and CTLA-4 were assessed. Interestingly, the expression level of CTLA-4 was significantly positively correlated with the GTscore, whereas that of GZMH and GZMK was significantly negatively correlated with the GTscore. CTLA-4 is a receptor that suppresses the T-cell function (Chen and Flies, 2013). Both GZMH and GZMK play a critical role in driving NK cell cytotoxicity and memory CD8 T-cell differentiation (Fellows et al., 2007; Harari et al., 2009). These results may provide a useful reference for us to further explore the immune microenvironment of NB.

As GD2-positive tumors are more responsive to anti-GD2 immunotherapy, it is necessary to screen GD2 molecular phenotypes on tumor cells. Recently, Sorokin et al. built a gene expression signature to predict the GD2-positive tumor phenotypes (Sorokin et al., 2020). Here, we investigated the ability of the GTscore to distinguish the GD2-positive phenotype in the four cases. Our analysis demonstrated that our GTscore displays good predictive ability for a putative GD2-positive phenotype. This model is helpful to predict the GD2 phenotype in clinical specimens from cancer patients. NSE can be useful for the diagnosis, prognosis, staging, and treatment of related neuroendocrine tumors (Faias et al., 2021). We found that our GTscore was positively correlated with NSE expression levels, indicating that the GTscore we developed could be used as a marker of neuroendocrine tumors.

Recently, immunotherapy has provided a promising strategy for NB (Kholodenko et al., 2018). However, only 20% to 30% of high-risk NB patients could benefit from immunotherapy (Voeller and Sondel, 2019). The current problem faced by researchers is how to identify patients who respond to immunotherapy. Although some studies have been done in this respect, there is still a lack of reliable markers. In the current study, we investigated the ability of the GTscore to predict the response of skin melanoma patients to anti–PD-1 immunotherapy. Our results indicated that the GTscore could reliably predict the anti–PD-1 therapy effect. In addition, the survival time of the low GTscore group was longer than that of the high GTscore group. On the basis of previous studies (Li et al., 1995), we hypothesize that cancer cells in the high GTscore subtype secrete the immunosuppressive molecule sGD2 (Ladisch et al., 1987), which inhibits T-cell proliferation and therefore likely contributes to immunotherapy resistance. However, the ability of the GTscore to predict the anti-GD2 immunotherapeutic response needs further verification.

There are several deficiencies in our research. First, another independent NB cohort that underwent GT assays and received anti-GD2 therapy with follow-up data will be required to validate our observation. Second, more laboratory experiments are needed to verify the detailed biological mechanism of the identified genes.

Conclusion

We constructed a GTscore, which could provide a favorable prognostic value in NB. Furthermore, we confirmed that the GTscore had a good predictive ability for putative GD2-positive molecular phenotype and NSE expression levels. Finally, the GTscore could accurately predict the efficacy of immunotherapy in melanoma. Nevertheless, large-scale, multicenter, and prospective studies are necessary to confirm our results.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/geo/, GSE49710, GSE90689, GSE142293 https://ocg.cancer.gov/programs/target, TARGET https://portal.gdc.cancer.gov/, TCGA.

Author Contributions

YS, LH, BS, and QZ: design, analysis and interpretation of data, drafting of the manuscript; LH: research direction and critical revision of the manuscript; QZ: research direction, critical revision of the manuscript, administrative support, obtaining funding, supervision. All authors read and agreed to this final manuscript.

Funding

This work was supported by grants from the National Key Research and Development Program of China (2018YFC1313000, 2018YFC1313001, and 2018YFC1313004).

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.

Acknowledgments

We appreciate the generosity of the studies from the TAGA, TARGET, and GEO groups for sharing large sets of data. We thank all the participants who were involved in this study and for providing editorial assistance.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.769580/full#supplementary-material

References

Aloia, A., Petrova, E., Tomiuk, S., Bissels, U., Déas, O., Saini, M., et al. (2015). The Sialyl-Glycolipid Stage-specific Embryonic Antigen 4 marks a Subpopulation of Chemotherapy-Resistant Breast Cancer Cells with Mesenchymal Features. Breast Cancer Res. 17 (1), 146. doi:10.1186/s13058-015-0652-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Asgharzadeh, S., Salo, J. A., Ji, L., Oberthuer, A., Fischer, M., Berthold, F., et al. (2012). Clinical Significance of Tumor-Associated Inflammatory Cells in Metastatic Neuroblastoma. Jco 30 (28), 3525–3532. doi:10.1200/JCO.2011.40.9169

CrossRef Full Text | Google Scholar

Bagchi, S., Yuan, R., and Engleman, E. G. (2021). Immune Checkpoint Inhibitors for the Treatment of Cancer: Clinical Impact and Mechanisms of Response and Resistance. Annu. Rev. Pathol. Mech. Dis. 16, 223–249. doi:10.1146/annurev-pathol-042020-042741

CrossRef Full Text | Google Scholar

Bakhoum, M. F., and Esmaeli, B. (2019). Molecular Characteristics of Uveal Melanoma: Insights from the Cancer Genome Atlas (TCGA) Project. Cancers 11 (8), 1061. doi:10.3390/cancers11081061

PubMed Abstract | CrossRef Full Text | Google Scholar

Balis, F. M., Busch, C. M., Desai, A. V., Hibbitts, E., Naranjo, A., Bagatell, R., et al. (2020). The Ganglioside G D2 as a Circulating Tumor Biomarker for Neuroblastoma. Pediatr. Blood Cancer 67 (1), e28031. doi:10.1002/pbc.28031

PubMed Abstract | CrossRef Full Text | Google Scholar

Berois, N., Gattolliat, C.-H., Barrios, E., Capandeguy, L., Douc-Rasy, S., Valteau-Couanet, D., et al. (2013). GALNT9 Gene Expression Is a Prognostic Marker in Neuroblastoma Patients. Clin. Chem. 59 (1), 225–233. doi:10.1373/clinchem.2012.192328

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhat, S. A., Mir, M. U. R., Majid, S., Hassan, T., Rehman, M. U., and Kuchy, S. (2018). Diagnostic Utility of Glycosyltransferase mRNA Expression in Gastric Cancer. Hematology/Oncology Stem Cel Ther. 11 (3), 158–168. doi:10.1016/j.hemonc.2018.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Bindea, G., Mlecnik, B., Tosolini, M., Kirilovsky, A., Waldner, M., Obenauf, A. C., et al. (2013). Spatiotemporal Dynamics of Intratumoral Immune Cells Reveal the Immune Landscape in Human Cancer. Immunity 39 (4), 782–795. doi:10.1016/j.immuni.2013.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Bocca, P., Di Carlo, E., Caruana, I., Emionite, L., Cilli, M., De Angelis, B., et al. (2017). Bevacizumab-mediated Tumor Vasculature Remodelling Improves Tumor Infiltration and Antitumor Efficacy of GD2-CAR T Cells in a Human Neuroblastoma Preclinical Model. Oncoimmunology 7 (1), e1378843. doi:10.1080/2162402X.2017.1378843

PubMed Abstract | CrossRef Full Text | Google Scholar

Braga, F., Ferraro, S., Mozzi, R., Dolci, A., and Panteghini, M. (2013). Biological Variation of Neuroendocrine Tumor Markers Chromogranin A and Neuron-specific Enolase. Clin. Biochem. 46 (1-2), 148–151. doi:10.1016/j.clinbiochem.2012.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Burchell, J. M., Beatson, R., Graham, R., Taylor-Papadimitriou, J., and Tajadura-Ortega, V. (2018). O-linked Mucin-type Glycosylation in Breast Cancer. Biochem. Soc. Trans. 46 (4), 779–788. doi:10.1042/BST20170483

PubMed Abstract | CrossRef Full Text | Google Scholar

Cavdarli, S., Groux-Degroote, S., and Delannoy, P. (2019). Gangliosides: The Double-Edge Sword of Neuro-Ectodermal Derived Tumors. Biomolecules 9 (8), 311. doi:10.3390/biom9080311

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, H.-H., Lu, M.-Y., Yang, Y.-L., Chou, S.-W., Lin, D.-T., Lin, K.-H., et al. (2020). The Prognostic Roles of and Correlation between ALK and MYCN Protein Expression in Neuroblastoma. J. Clin. Pathol. 73 (3), 154–161. doi:10.1136/jclinpath-2019-206063

CrossRef Full Text | Google Scholar

Chatterjee, S. B., Hou, J., Bandaru, V. V. R., Pezhouh, M. K., Syed Rifat Mannan, A. A., and Sharma, R. (2019). Lactosylceramide Synthase β-1,4-GalT-V: A Novel Target for the Diagnosis and Therapy of Human Colorectal Cancer. Biochem. Biophysical Res. Commun. 508 (2), 380–386. doi:10.1016/j.bbrc.2018.11.149

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, L., and Flies, D. B. (2013). Molecular Mechanisms of T Cell Co-stimulation and Co-inhibition. Nat. Rev. Immunol. 13 (4), 227–242. doi:10.1038/nri3405

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheresh, D. A., Rosenberg, J., Mujoo, K., Hirschowitz, L., and Reisfeld, R. A. (1986). Biosynthesis and Expression of the Disialoganglioside GD2, a Relevant Target Antigen on Small Cell Lung Carcinoma for Monoclonal Antibody-Mediated Cytolysis. Cancer Res. 46 (10), 5112–5118.

PubMed Abstract | Google Scholar

Ding, W., Xu, X., Qian, Y., Xue, W., Wang, Y., Du, J., et al. (2018). Prognostic Value of Tumor-Infiltrating Lymphocytes in Hepatocellular Carcinoma. Medicine (Baltimore) 97 (50), e13301. doi:10.1097/MD.0000000000013301

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, S., Wang, Z., Huang, B., Zhang, J., Ge, Y., Fan, Q., et al. (2017). Bioinformatics Insight into Glycosyltransferase Gene Expression in Gastric Cancer: POFUT1 Is a Potential Biomarker. Biochem. Biophysical Res. Commun. 483 (1), 171–177. doi:10.1016/j.bbrc.2016.12.172

CrossRef Full Text | Google Scholar

Faias, S., Prazeres, S., Cunha, M., Pereira, L., Roque, R., Chaves, P., et al. (2020). Chromogranin A and NSE in Cystic Pancreatic Neuroendocrine Tumors. Clin. Res. Hepatol. Gastroenterol. 45 (4), 101601. doi:10.1016/j.clinre.2020.101601

PubMed Abstract | CrossRef Full Text | Google Scholar

Fardin, P., Barla, A., Mosci, S., Rosasco, L., Verri, A., Versteeg, R., et al. (2010). A Biology-Driven Approach Identifies the Hypoxia Gene Signature as a Predictor of the Outcome of Neuroblastoma Patients. Mol. Cancer 9, 185. doi:10.1186/1476-4598-9-185

PubMed Abstract | CrossRef Full Text | Google Scholar

Fellows, E., Gil-Parrado, S., Jenne, D. E., and Kurschus, F. C. (2007). Natural Killer Cell-Derived Human Granzyme H Induces an Alternative, Caspase-independent Cell-Death Program. Blood 110 (2), 544–552. doi:10.1182/blood-2006-10-051649

PubMed Abstract | CrossRef Full Text | Google Scholar

Fisher, J. P. H., Flutter, B., Wesemann, F., Frosch, J., Rossig, C., Gustafsson, K., et al. (2016). Effective Combination Treatment of GD2-Expressing Neuroblastoma and Ewing's Sarcoma Using Anti-GD2 ch14.18/CHO Antibody with Vγ9Vδ2+ γδT Cells. Oncoimmunology 5 (1), e1025194. doi:10.1080/2162402X.2015.1025194

PubMed Abstract | CrossRef Full Text | Google Scholar

Galvan, M., Tsuboi, S., Fukuda, M., and Baum, L. G. (2000). Expression of a Specific Glycosyltransferase Enzyme Regulates T Cell Death Mediated by Galectin-1. J. Biol. Chem. 275 (22), 16730–16737. doi:10.1074/jbc.M001117200

CrossRef Full Text | Google Scholar

Giovanella, L., Ceriani, L., Bandera, M., Beghe, B., and Roncari, G. (1995). Evaluation of the Serum Markers CEA, NSE, TPS and CYFRA 21.1 in Lung Cancer. Int. J. Biol. Markers 10 (3), 156–160. doi:10.1177/172460089501000305

PubMed Abstract | CrossRef Full Text | Google Scholar

Groux-Degroote, S., Guérardel, Y., and Delannoy, P. (2017). Gangliosides: Structures, Biosynthesis, Analysis, and Roles in Cancer. Chembiochem 18 (13), 1146–1154. doi:10.1002/cbic.201600705

PubMed Abstract | CrossRef Full Text | Google Scholar

Gupta, R., Leon, F., Thompson, C. M., Nimmakayala, R., Karmakar, S., Nallasamy, P., et al. (2020). Global Analysis of Human Glycosyltransferases Reveals Novel Targets for Pancreatic Cancer Pathogenesis. Br. J. Cancer 122 (11), 1661–1672. doi:10.1038/s41416-020-0772-3

CrossRef Full Text | Google Scholar

Hamilton, W. B., Helling, F., Lloyd, K. O., and Livingston, P. O. (1993). Ganglioside Expression on Human Malignant Melanoma Assessed by Quantitative Immune Thin-Layer Chromatography. Int. J. Cancer 53 (4), 566–573. doi:10.1002/ijc.2910530407

CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Harari, A., Enders, F. B., Cellerai, C., Bart, P.-A., and Pantaleo, G. (2009). Distinct Profiles of Cytotoxic Granules in Memory CD8 T Cells Correlate with Function, Differentiation Stage, and Antigen Exposure. J. Virol. 83 (7), 2862–2871. doi:10.1128/JVI.02528-08

CrossRef Full Text | Google Scholar

Hu, H., Li, Z., Zhou, Y., Zhang, Y., Zhao, L., Zhao, W., et al. (2019). GLT8D1 Overexpression as a Novel Prognostic Biomarker in Human Cutaneous Melanoma. Melanoma Res. 29 (6), 612–620. doi:10.1097/CMR.0000000000000631

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, Q., Miller, M. R., Schappet, J., and Henry, M. D. (2015). The Glycosyltransferase LARGE2 Is Repressed by Snail and ZEB1 in Prostate Cancer. Cancer Biol. Ther. 16 (1), 125–136. doi:10.4161/15384047.2014.987078

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, S., Cai, N., Pacheco, P. P., Narrandes, S., Wang, Y., and Xu, W. (2018). Applications of Support Vector Machine (SVM) Learning in Cancer Genomics. Cgp 15 (1), 41–51. doi:10.21873/cgp.20063

PubMed Abstract | CrossRef Full Text | Google Scholar

Iwasawa, T., Zhang, P., Ohkawa, Y., Momota, H., Wakabayashi, T., Ohmi, Y., et al. (2018). Enhancement of Malignant Properties of Human Glioma Cells by Ganglioside GD3/GD2. Int. J. Oncol. 52 (4), 1255–1266. doi:10.3892/ijo.2018.4266

CrossRef Full Text | Google Scholar

Jiang, J., and Gu, J. (2010). β1,4-Galactosyltransferase V. Methods Enzymol. 479, 3–23. doi:10.1016/S0076-6879(10)79001-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, Y., and Yang, Y. (2019). Bioinformatics-based Discovery of PYGM and TNNC2 as Potential Biomarkers of Head and Neck Squamous Cell Carcinoma. Biosci. Rep. 39 (7). doi:10.1042/BSR20191612

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, B. D., Jing, W., and Orentas, R. J. (2007). CD25+ Regulatory T Cell Inhibition Enhances Vaccine-Induced Immunity to Neuroblastoma. J. Immunother. 30 (2), 203–214. doi:10.1097/01.cji.0000211336.91513.dd

CrossRef Full Text | Google Scholar

Joshi, S. (2020). Targeting the Tumor Microenvironment in Neuroblastoma: Recent Advances and Future Directions. Cancers 12 (8), 2057. doi:10.3390/cancers12082057

PubMed Abstract | CrossRef Full Text | Google Scholar

Kholodenko, I. V., Kalinovsky, D. V., Doronin, , Deyev, S. M., and Kholodenko, R. V. (2018). Neuroblastoma Origin and Therapeutic Targets for Immunotherapy. J. Immunol. Res. 2018, 1–25. doi:10.1155/2018/7394268

PubMed Abstract | CrossRef Full Text | Google Scholar

Ladisch, S., Becker, H., and Ulsh, L. (1992). Immunosuppression by Human Gangliosides: I. Relationship of Carbohydrate Structure to the Inhibition of T Cell Responses. Biochim. Biophys. Acta (Bba) - Lipids Lipid Metab. 1125 (2), 180–188. doi:10.1016/0005-2760(92)90043-u

CrossRef Full Text | Google Scholar

Li, R., Villacreses, N., and Ladisch, S. (1995). Human Tumor Gangliosides Inhibit Murine Immune Responses In Vivo. Cancer Res. 55 (2), 211–214.

PubMed Abstract | Google Scholar

Lin, B., Du, L., Li, H., Zhu, X., Cui, L., and Li, X. (2020). Tumor-infiltrating Lymphocytes: Warriors Fight against Tumors Powerfully. Biomed. Pharmacother. 132, 110873. doi:10.1016/j.biopha.2020.110873

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, K. X., and Joshi, S. (2020). "Re-educating" Tumor Associated Macrophages as a Novel Immunotherapy Strategy for Neuroblastoma. Front. Immunol. 11, 1947. doi:10.3389/fimmu.2020.01947

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X., Zhang, Y., Han, Y., Lu, W., Yang, J., Tian, J., et al. (2020). Overexpression of GLT1D1 Induces Immunosuppression through Glycosylation of PD‐L1 and Predicts Poor Prognosis in B‐cell Lymphoma. Mol. Oncol. 14 (5), 1028–1044. doi:10.1002/1878-0261.12664

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Wondimu, A., Yan, S., Bobb, D., and Ladisch, S. (2014). Tumor Gangliosides Accelerate Murine Tumor Angiogenesis. Angiogenesis 17 (3), 563–571. doi:10.1007/s10456-013-9403-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Luen, S. J., Savas, P., Fox, S. B., Salgado, R., and Loi, S. (2017). Tumour-infiltrating Lymphocytes and the Emerging Role of Immunotherapy in Breast Cancer. Pathology 49 (2), 141–155. doi:10.1016/j.pathol.2016.10.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Macheri, B., Lee, W. F., and Westrick, M. (1982). Glycosphingolipids of normal and Leukemic Human Leukocytes. Mol. Cel Biochem 47 (2), 81–95. doi:10.1007/BF00234409

PubMed Abstract | CrossRef Full Text | Google Scholar

Mansoori, M., Roudi, R., Abbasi, A., Abolhasani, M., Abdi Rad, I., Shariftabrizi, A., et al. (2019). High GD2 Expression Defines Breast Cancer Cells with Enhanced Invasiveness. Exp. Mol. Pathol. 109, 25–35. doi:10.1016/j.yexmp.2019.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Maris, J. M., Hogarty, M. D., Bagatell, R., and Cohn, S. L. (2007). Neuroblastoma. The Lancet 369 (9579), 2106–2120. doi:10.1016/S0140-6736(07)60983-0

CrossRef Full Text | Google Scholar

Mehta, K. A., Patel, K. A., Pandya, S. J., and Patel, P. S. (2020). "Aberrant Sialylation Plays a Significant Role in Oral Squamous Cell Carcinoma Progression". J. Oral Pathol. Med. 49 (3), 253–259. doi:10.1111/jop.12976

CrossRef Full Text | Google Scholar

Mennel, H. D., Bosslet, K., Geissel, H., and Bauer, B. L. (2000). Immunohistochemically Visualized Localisation of Gangliosides Glac2 (GD3) and Gtri2 (GD2) in Cells of Human Intracranial Tumors. Exp. Toxicologic Pathol. 52 (4), 277–285. doi:10.1016/S0940-2993(00)80046-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Metelitsa, L. S., Wu, H.-W., Wang, H., Yang, Y., Warsi, Z., Asgharzadeh, S., et al. (2004). Natural Killer T Cells Infiltrate Neuroblastomas Expressing the Chemokine CCL2. J. Exp. Med. 199 (9), 1213–1221. doi:10.1084/jem.20031462

CrossRef Full Text | Google Scholar

Mina, M., Boldrini, R., Citti, A., Romania, P., D'Alicandro, V., De Ioris, M., et al. (2015). Tumor-infiltrating T Lymphocytes Improve Clinical Outcome of Therapy-Resistant Neuroblastoma. Oncoimmunology 4 (9), e1019981. doi:10.1080/2162402X.2015.1019981

PubMed Abstract | CrossRef Full Text | Google Scholar

Mujoo, K., Cheresh, D. A., Yang, H. M., and Reisfeld, R. A. (1987). Disialoganglioside GD2 on Human Neuroblastoma Cells: Target Antigen for Monoclonal Antibody-Mediated Cytolysis and Suppression of Tumor Growth. Cancer Res. 47 (4), 1098–1104.

PubMed Abstract | Google Scholar

Narvaez, C. J., Grebenc, D., Balinth, S., and Welsh, J. E. (2020). Vitamin D Regulation of HAS2, Hyaluronan Synthesis and Metabolism in Triple Negative Breast Cancer Cells. J. Steroid Biochem. Mol. Biol. 201, 105688. doi:10.1016/j.jsbmb.2020.105688

CrossRef Full Text | Google Scholar

Ohmi, Y., Kambe, M., Ohkawa, Y., Hamamura, K., Tajima, O., Takeuchi, R., et al. (2018). Differential Roles of Gangliosides in Malignant Properties of Melanomas. PLoS One 13 (11), e0206881. doi:10.1371/journal.pone.0206881

PubMed Abstract | CrossRef Full Text | Google Scholar

Orsi, G., Barbolini, M., Ficarra, G., Tazzioli, G., Manni, P., Petrachi, T., et al. (2017). GD2 Expression in Breast Cancer. Oncotarget 8 (19), 31592–31600. doi:10.18632/oncotarget.16363

PubMed Abstract | CrossRef Full Text | Google Scholar

Pérel, Y., Valteau-Couanet, D., Michon, J., Lavrand, F., Coze, C., Bergeron, C., et al. (2004). Le pronostic du neuroblastome de l'enfant. Méthodes d'étude et utilisation en pratique clinique. Arch. de Pédiatrie 11 (7), 834–842. doi:10.1016/j.arcped.2004.02.022

CrossRef Full Text | Google Scholar

Pérez-Garay, M., Arteta, B., Pagès, L., de Llorens, R., de Bolòs, C., Vidal-Vanaclocha, F., et al. (2010). α2,3-Sialyltransferase ST3Gal III Modulates Pancreatic Cancer Cell Motility and Adhesion In Vitro and Enhances its Metastatic Potential In Vivo. PLoS One 5 (9), e12524. doi:10.1371/journal.pone.0012524

PubMed Abstract | CrossRef Full Text | Google Scholar

Pick, J., Anh-Tuan, N., Mód, A., and Hollán, S. (1986). Ganglioside Composition in Common Acute Lymphoblastic Leukaemia. Haematologia (Budap) 19 (1), 33–9.

PubMed Abstract | Google Scholar

Poon, V. I., Roth, M., Piperdi, S., Geller, D., Gill, J., Rudzinski, E. R., et al. (2015). Ganglioside GD2 Expression Is Maintained upon Recurrence in Patients with Osteosarcoma. Clin. Sarcoma Res. 5 (1), 4. doi:10.1186/s13569-014-0020-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Potapenko, I. O., Haakensen, V. D., Lüders, T., Helland, Å., Bukholm, I., Sørlie, T., et al. (2010). Glycan Gene Expression Signatures in normal and Malignant Breast Tissue; Possible Role in Diagnosis and Progression. Mol. Oncol. 4 (2), 98–118. doi:10.1016/j.molonc.2009.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Rakhshandehroo, M., van Eijkeren, R. J., Gabriel, T. L., de Haar, C., Gijzel, S. M. W., Hamers, N., et al. (2019). Adipocytes Harbor a Glucosylceramide Biosynthesis Pathway Involved in iNKT Cell Activation. Biochim. Biophys. Acta (Bba) - Mol. Cel Biol. Lipids 1864 (8), 1157–1167. doi:10.1016/j.bbalip.2019.04.016

CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Rovigatti, U. (2021). The Glycosphingolipid GD2 as an Effective but Enigmatic Target of Passive Immunotherapy in Children with Aggressive Neuroblastoma (HR-NBL). Cancer Lett. 503, 220–230. doi:10.1016/j.canlet.2020.11.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Saito, S., Aoki, H., Ito, A., Ueno, S., Wada, T., Mitsuzuka, K., et al. (2003). Human α2,3-Sialyltransferase (ST3Gal II) Is a Stage-specific Embryonic Antigen-4 Synthase. J. Biol. Chem. 278 (29), 26474–26479. doi:10.1074/jbc.M213223200

PubMed Abstract | CrossRef Full Text | Google Scholar

Schultz, M. J., Holdbrooks, A. T., Chakraborty, A., Grizzle, W. E., Landen, C. N., Buchsbaum, D. J., et al. (2016). The Tumor-Associated Glycosyltransferase ST6Gal-I Regulates Stem Cell Transcription Factors and Confers a Cancer Stem Cell Phenotype. Cancer Res. 76 (13), 3978–3988. doi:10.1158/0008-5472.CAN-15-2834

PubMed Abstract | CrossRef Full Text | Google Scholar

Schulz, G., Cheresh, D. A., Varki, N. M., Yu, A., Staffileno, L. K., and Reisfeld, R. A. (1984). Detection of Ganglioside GD2 in Tumor Tissues and Sera of Neuroblastoma Patients. Cancer Res. 44 (12 Pt 1), 5914–5920.

PubMed Abstract | Google Scholar

Seko, A., Kataoka, F., Aoki, D., Sakamoto, M., Nakamura, T., Hatae, M., et al. (2009). β1,3-Galactosyltransferases-4/5 Are Novel Tumor Markers for Gynecological Cancers. Tumor Biol. 30 (1), 43–50. doi:10.1159/000203129

CrossRef Full Text | Google Scholar

Shine, B. S., Hungerford, J., Vaghela, B., and Sheraidah, G. A. (1990). Electrophoretic Assessment of Aqueous and Serum Neurone-specific Enolase in Retinoblastoma and Ocular Malignant Melanoma. Br. J. Ophthalmol. 74 (7), 427–430. doi:10.1136/bjo.74.7.427

CrossRef Full Text | Google Scholar

Sokol, E., Desai, A. V., Applebaum, M. A., Valteau-Couanet, D., Park, J. R., Pearson, A. D. J., et al. (2020). Age, Diagnostic Category, Tumor Grade, and Mitosis-Karyorrhexis Index Are Independently Prognostic in Neuroblastoma: An INRG Project. Jco 38 (17), 1906–1918. doi:10.1200/JCO.19.03285

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, J., Zhao, Q., Xu, Y., and Zhu, L. (2020). A Signature of 29 Immune-Related Genes Pairs to Predict Prognosis in Patients with Neuroblastoma. Int. Immunopharmacology 88, 106994. doi:10.1016/j.intimp.2020.106994

CrossRef Full Text | Google Scholar

Sorice, M., Longo, A., Garofalo, T., Mattei, V., Misasi, R., and Pavan, A. (2003). Role of GM3-Enriched Microdomains in Signal Transduction Regulation in T Lymphocytes. Glycoconj J. 20 (1), 63–70. doi:10.1023/B:GLYC.0000018018.29488.c6

CrossRef Full Text | Google Scholar

Sorokin, M., Kholodenko, I., Kalinovsky, D., Shamanskaya, T., Doronin, I., Konovalov, D., et al. (2020). RNA Sequencing-Based Identification of Ganglioside GD2-Positive Cancer Phenotype. Biomedicines 8 (6), 142. doi:10.3390/biomedicines8060142

PubMed Abstract | CrossRef Full Text | Google Scholar

Suzuki, M., and Cheung, N.-K. V. (2015). Disialoganglioside GD2 as a Therapeutic Target for Human Diseases. Expert Opin. Ther. Targets 19 (3), 349–362. doi:10.1517/14728222.2014.986459

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Franceschini, A., Wyder, S., Forslund, K., Heller, D., Huerta-Cepas, J., et al. (2015). STRING V10: Protein-Protein Interaction Networks, Integrated over the Tree of Life. Nucleic Acids Res. 43 (Database issue), D447–D452. doi:10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, L., Chen, X., Zhang, X., Guo, Y., Su, J., Zhang, J., et al. (2019). N-glycosylation in Progression of Skin Cancer. Med. Oncol. 36 (6), 50. doi:10.1007/s12032-019-1270-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, W., Li, M., Qi, X., and Li, J. (2020). β1,4-Galactosyltransferase V Modulates Breast Cancer Stem Cells through Wnt/β-Catenin Signaling Pathway. Cancer Res. Treat. 52 (4), 1084–1102. doi:10.4143/crt.2020.093

PubMed Abstract | CrossRef Full Text | Google Scholar

Terzic, T., Cordeau, M., Herblot, S., Teira, P., Cournoyer, S., Beaunoyer, M., et al. (2018). Expression of Disialoganglioside (GD2) in Neuroblastic Tumors: A Prognostic Value for Patients Treated with Anti-GD2 Immunotherapy. Pediatr. Dev. Pathol. 21 (4), 355–362. doi:10.1177/1093526617723972

PubMed Abstract | CrossRef Full Text | Google Scholar

Tibshirani, R. (1997). The Lasso Method for Variable Selection in the Cox Model. Statist. Med. 16 (4), 385–395. doi:10.1002/(sici)1097-0258(19970228)16:4<385:aid-sim380>3.0.co

CrossRef Full Text | Google Scholar

Troschke-Meurer, S., Siebert, N., Marx, M., Zumpe, M., Ehlert, K., Mutschlechner, O., et al. (2019). Low CD4⁺/CD25⁺/CD127⁻ Regulatory T Cell- and High INF-γ Levels Are Associated with Improved Survival of Neuroblastoma Patients Treated with Long-Term Infusion of ch14.18/CHO Combined with Interleukin-2. Oncoimmunology 8 (12), 1661194. doi:10.1080/2162402X.2019.1661194

PubMed Abstract | CrossRef Full Text | Google Scholar

Vanichapol, T., Chutipongtanate, S., Anurathapan, U., and Hongeng, S. (2018). Immune Escape Mechanisms and Future Prospects for Immunotherapy in Neuroblastoma. Biomed. Res. Int. 2018, 1–11. doi:10.1155/2018/1812535

PubMed Abstract | CrossRef Full Text | Google Scholar

Voeller, J., and Sondel, P. M. (2019). Advances in Anti-GD2 Immunotherapy for Treatment of High-Risk Neuroblastoma. J. Pediatr. Hematol. Oncol. 41 (3), 163–169. doi:10.1097/MPH.0000000000001369

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Wang, X., Xu, L., Zhang, J., and Cao, H. (2020). Prognostic Significance of MYCN Related Genes in Pediatric Neuroblastoma: a Study Based on TARGET and GEO Datasets. BMC Pediatr. 20 (1), 314. doi:10.1186/s12887-020-02219-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Li, S., Yu, X., Han, Y., Wu, Y., Wang, S., et al. (2019). α2,6-Sialylation Promotes Immune Escape in Hepatocarcinoma Cells by Regulating T Cell Functions and CD147/MMP Signaling. J. Physiol. Biochem. 75 (2), 199–207. doi:10.1007/s13105-019-00674-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Whiteside, T. L. (2008). The Tumor Microenvironment and its Role in Promoting Tumor Growth. Oncogene 27 (45), 5904–5912. doi:10.1038/onc.2008.271

PubMed Abstract | CrossRef Full Text | Google Scholar

Wienke, J., Dierselhuis, M. P., Tytgat, G. A. M., Künkele, A., Nierkens, S., and Molenaar, J. J. (2021). The Immune Landscape of Neuroblastoma: Challenges and Opportunities for Novel Therapeutic Strategies in Pediatric Oncology. Eur. J. Cancer 144, 123–150. doi:10.1016/j.ejca.2020.11.014

CrossRef Full Text | Google Scholar

Wu, Z.-L., Ladisch, S., Feig, S., Ulsh, L., Schwartz, E., Floutsis, G., et al. (1987). Shedding of GD2 Ganglioside by Human Neuroblastoma. Int. J. Cancer 39 (1), 73–76. doi:10.1002/ijc.2910390113

CrossRef Full Text | Google Scholar

Wu, Z. L., Schwartz, E., Seeger, R., and Ladisch, S. (1986). Expression of GD2 Ganglioside by Untreated Primary Human Neuroblastomas. Cancer Res. 46 (1), 440–443.

PubMed Abstract | Google Scholar

Xu, S., Zhang, S., Chen, C., Yan, J., Cai, M., Zhu, X., et al. (2002). Over-expression of Beta-1,4-Galactosyltransferase V Increases the Growth of Astrocytoma Cell Line. J. Exp. Clin. Cancer Res. 21 (3), 409–414.

Google Scholar

Yoshida, S., Fukumoto, S., Kawaguchi, H., Sato, S., Ueda, R., and Furukawa, K. (2001). Ganglioside G(D2) in Small Cell Lung Cancer Cell Lines: Enhancement of Cell Proliferation and Mediation of Apoptosis. Cancer Res. 61 (10), 4244–4252.

PubMed Abstract | Google Scholar

Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16 (5), 284–287. doi:10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

Zeltzer, P., Parma, A., Dalton, A., Siegel, S., Marangos, P., Sather, H., et al. (1983). Raised Neuron-specific Enolase in Serum of Children with Metastatic Neuroblastoma. The Lancet 322 (8346), 361–363. doi:10.1016/s0140-6736(83)90342-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, T., Wang, F., Wu, J.-Y., Qiu, Z.-C., Wang, Y., Liu, F., et al. (2018). Clinical Correlation of B7-H3 and B3GALT4 with the Prognosis of Colorectal Cancer. Wjg 24 (31), 3538–3546. doi:10.3748/wjg.v24.i31.3538

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, W., Yu, Y., Hertwig, F., Thierry-Mieg, J., Zhang, W., Thierry-Mieg, D., et al. (2015). Comparison of RNA-Seq and Microarray-Based Models for Clinical Endpoint Prediction. Genome Biol. 16, 133. doi:10.1186/s13059-015-0694-1

PubMed Abstract | CrossRef Full Text | Google Scholar