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.

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.

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.

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).

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: GTscore  =  β 1 ∗X 1 + β 2 ∗X 2 + … βi∗Xi  … +  β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.

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.

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).

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.

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/).

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.

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.

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.

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).

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).

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.

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).

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.

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.