Exploring the role of neutrophil extracellular traps in neuroblastoma: identification of molecular subtypes and prognostic implications

Abstract

Background: Cancer cells induce neutrophil extracellular traps (NETs) to promote tumor progression and metastasis. However, only a few studies have focused on the role of NETs in Neuroblastoma (NB).

Methods: First, based on the expression of NET-related genes, consensus clustering analysis was conducted to cluster NB samples into different subtypes. Differential analysis was performed to identify DEGs between subtypes. Functional items and related pathways of DEGs were identified using enrichment analysis. Univariate Cox analysis and the LASSO algorithm were used to identify biomarkers for prognosis. Furthermore, independent prognostic analysis was performed. Immune infiltration analysis was performed to identify differential immune cells. Finally, the verification of prognostic model genes were taken by the immunohistochemical staining and quantitative real-time PCR.

Results: Consensus clustering analysis demonstrated that NB samples were clustered into two subtypes. There were 125 DEGs between the two subtypes of NB. Moreover, the enrichment analysis results showed that the DEGs were mainly associated with 'external side of plasma membrane,' 'immune receptor activity' 'regulation of leukocyte migration' GO items. There were also several GO items related to neutrophils, such as regulation of neutrophil migration and differentiation. KEGG pathways revealed that the DEGs were correlated with in immunity-related activities, including 'Complement and coagulation cascades,' 'Neutrophil extracellular trap formation, 'T cell receptor signaling pathway,' 'PD-L1 expression and PD-1 checkpoint pathway in cancer' and so on. A total of five biomarkers,[Selenoprotein P1 (SEPP1), Fibrinogen-like protein 2 (FGL2), NK cell lectin-like receptor K1 (KLRK1), ATP-binding cassette transporters 6(ABCA6) and Galectins(GAL)], were screened, and a risk model based on the biomarkers was created. Furthermore, a nomogram for forecasting the survival rates of patients with NB was established based on the risk score, age at diagnosis, and MYCN status. Eight differential immune cells (CD8 + T cells, resting mast cells, etc.) were acquired between the two risk subgroups. The expression levels of five prognostic model genes at the protein and mRNA were verified and all results were consistent with the results of our bioinformatics analysis.

Conclusion: We initially found that five NET-related genes were significantly differentially expressed in NETs-associated molecular isoforms and two Netrg molecular isoforms were found to be associated with poorer prognosis. This stratification might provide insight into the prediction of prognosis and ideal immunotherapy strategies for patients with NB. However, we also noted that the formation of NETs is a complex biological process involving the regulation of multiple cytokines and cellular interactions. Therefore, the exact roles of these genes and their specific mechanisms in the formation of NETs and the development of NB still need to be further investigated.

Keywords: immune; inflammation; neuroblastoma; neutrophil extracellular traps; prognosis.

Figure 1

Identification of NETs -related subtypes (A) Consensus matrix heat map based on the expression of 136 NETs-related genes using R packege ConsensusClusterPlus on the NBL samples in the training set. (B) Cumulative distribution of CDF to determine the optimal number of subtypes (C) Delta area plot showing the relative change of areas under the CDF curve (D) Principal component analysis graphs generated based on cluster 1 and cluster 2 in the training set after classification into 2 subtypes.

Figure 2

Identification of differentially expressed genes associated with NETs, were plotted using R package ggpubr and pheatmap, enrichment analysis of DEG, respectively. (A) The volcano maps, there were 125 differentially expressed genes between cluster 1 vs cluster 2 groups, of which 123 were up-regulated expressed genes and 2 were down-regulated expressed genes. (B) The heatmaps,the red dots on the graph indicate up-regulated genes, the blue dots indicate down-regulated genes, the black ones are non-significant, and the more skewed towards the upper left and right corners of the graph, the greater the multiplicity of differences in significance. (C) GO and (D) KEGG enrichment analyses,125 NETs-related differential genes for GO functional annotation.

Figure 3

Development of the NETs -related prognostic signature in the training set (A) Forest plot of 26 NETs-related genes correlated with OS. (B1-2) The LASSO Cox regression analysis was performed depending on the optimal λ value. (C) Distribution of prognostic index in different risk groups, survival status of patients in different risk groups and the heatmap of prognostic gene. (D) The Kaplan-Meier curves showed that significant differences were identified for OS between these two risk groups. (E) ROC curves evaluating the sensitivity and specificity of the NETs-based prognostic model.

Figure 4

Validation of the NETs -based signature in the testing set.(A) Distribution of TARGET prognostic index in high risk group (n = 108) and low risk group (n = 42), survival status of patients in high risk group (n = 108) and low risk group (n = 42) and the heatmap of prognostic gene. (B) The Kaplan-Meier curves showed that significant differences were identified for OS between high risk group (n = 108) and low risk group (n = 42). (C) ROC curves evaluating the sensitivity and specificity of the NETs-based prognostic model(TARGET). (D–G) The association between clinical traits and risk score of NB.

Figure 5

Independent prognostic analysis of risk scores and clinical parameters. (A, B) Forest plot of univariable and multivariate Cox model for survival in NB. “Age_at_diagnosis”: age at diagnosis; “MYCN_status”: the status of the MYCN gene in the tumour cells of a patient with neuroblastoma, a proto-oncogene associated with tumour growth and spread that is amplified in some neuroblastomas and is associated with aggressiveness and poor prognosis of the disease. “PROGRESSION”: in medicine and biology, the development or worsening of a disease or condition. “INSS_stage”: the abbreviation for international neuroblastoma staging system (INSS). INSS stages are as follows (St: Stage), Stage 1: The tumour is confined to the primary site, can be completely resected, and the patient is 18 months of age or older or the tumour diameter is less than or equal to 5 cm; Stage 2A: The tumour is partially resected and the ipsilateral residual tumour does not cross the midline; Stage 2B: partial resection of the tumour with ipsilateral residual tumour crossing the midline; Stage 3: Incomplete resection of the tumour due to invasion of major vascular structures or tumour crossing the midline with residual tumour; Stage 4: Tumour has distant metastases to other parts of the body; Stage 4S: Low risk of metastasis under certain conditions, usually seen in infants. (C) A nomogram to predict 1-year, 2-years, and 3-years survival rates in NB patients. (D) Nomogram-Predicted Probability of 1-3 year OS. (E1–3) Survival-dependent receiver operating characteristic curves for risk score, nomogram, and clinical pathological characteristics.

Figure 6

Immune-related characteristics in the low- and high-risk score groups. (A) Violin plots of stromal score, IMMUNE score, and ESTIMATE score for the training set samples in the high (n=51) and low (n=218) risk groups based on the training set using R package estimate. (B) The correlation between risk scores and stomal score, immune score, ESTIMATE score. (C) The distribution proportion of immune cell abundance in 269 sample. (D) Differences in immune cell infiltration between the low (n=218) and high (n=51) risk score groups. (E) The correlations of eight immune cells and risk scores. *p<0.05, ** p<0.01, *** p<0.001, **** p<0.0001, ns p>0.05.

Figure 7

Immunohistochemical results and column diagram of relative IOD values of five prognostic model genes (GAL, SEPP1, FGL2, KLRK1, and ABCA6) in the High-Risk (n=7) NB and Non High-Risk (n=25) NB samples.

Figure 8

Relative mRNA expressing of five prognostic model genes (GAL, SEPP1, FGL2, KLRK1, and ABCA6) expression in the High-Risk NB and Non High-Risk NB samples by the quantitative real-time PCR.