Tertiary lymphoid structures-driven immune infiltration patterns and their association with survival in neuroblastoma

Abstract

Background: Neuroblastoma (NB), a diverse childhood cancer, needs better prognostic markers for personalized treatment. The current clinical risk stratification system does not fully explain the high heterogeneity of tumor patients. Tertiary lymphoid structures (TLS), key in tumor immunity, may serve as new biomarkers, but their impact on NB prognosis is unclear.

Methods: We combined transcriptome data from NB cohorts GSE49710 and GSE62564, analyzing 37 TLS-related genes. A prognostic signature (CMLS) was created using machine learning and validated with Kaplan-Meier and receiver operating characteristic (ROC) curves. We also studied immune infiltration and gene expression patterns in NB tissues using single-cell sequencing and quantitative real-time polymerase chain reaction (qRT-PCR).

Results: A 6-gene TLS signature predicted better survival in NB patients. High levels of CCL2, CCL4, CCL21, CD200, CXCR3, and IGSF6 correlated with improved survival. The low-TLS risk group showed better event-free and overall survival. Immune analysis indicated a higher immune cell presence, especially cytotoxic T cells, in this group. Single-cell sequencing revealed lower TLS gene expression in refractory recurrence samples. CD200 downregulation reduced NB cell invasiveness and migration.

Conclusion: Our study demonstrates that TLS-related genes play a crucial role in NB prognosis, with a 6-gene TLS signature (CCL2, CCL4, CCL21, CD200, CXCR3, and IGSF6) serving as a promising prognostic biomarker for NB. CD200 may be a potential target for inhibiting the biological behavior of NB cells.

Keywords: Consensus machine learning-driven signature; Neuroblastoma; Single-cell sequencing; Tertiary lymphoid structures; Tumor microenvironment.

Figure 1. Screening of prognostic genes related to the tertiary lymphoid structure of neuroblastoma.

(A) Forest plot of univariate COX regression analysis based on the GSE49710 dataset. (B) Forest plot of univariate COX regression analysis based on the GSE62564 dataset.

Figure 2. Construction of consensus machine learning-based prognostic signature.

(A) A combination algorithm of 10 machine learning methods, the blue column represents the results for GSE49710, the green column corresponds to the results for GSE62564, and the third column corresponds to the C-index. (B) It displays the number of features each algorithm combination suggests to include. (C) For the constructed characteristic curve, the horizontal coordinate is the characteristic number and the vertical coordinate is the AUC value. (D) BP, CC and MF represent biological processes, cellular components and molecular functions, respectively.

Figure 3. Survival analysis of candidate genes in NB samples.

(A–F) was K-M analysis of six candidate genes in the GSE49710 dataset, respectively. (G–L) was the K-M analysis of six candidate genes in the GSE49710 dataset, respectively. The red curve represents the group with high gene expression and the blue curve represents the group with low gene expression.

Figure 4. Comprehensive analysis of the performance of TLS scoring model.

(A and B), (G and H) are K-M analysis curves of OS and EFS of TLS scoring model, respectively. The red curve represents NB groups with high ratings and the blue represents NB groups with low ratings. (C and D) and (I and J) are the diagnostic dependent ROC curves of the model respectively, and the gray shadow area corresponds to the AUC value. (E and F) and (K and L) were the time-dependent ROC curves of 0.5 and 1 year of the scoring model, respectively. Blue corresponds to 0.5 years, red corresponds to 1 year of survival.

Figure 5. Analysis of immune infiltration.

(A, B–D) and (E–G) correspond to the analysis results of TIMER, ESTIMATE and CIBERSORT, respectively. The red column represents the low rating group, and the blue represents the high rating group in (A–D). In (E), red represents the high rating group and green represents the low rating group. The size of the circle in the F diagram represents the degree of correlation. *, **, *** represent P < 0.05, P < 0.01, P < 0.001, respectively.

Figure 6. The correlation analysis between the TLS model and TME, hypoxia scores, vascular scores and tumor stemness.

(A–D) and (E–H) showed the relationship between TLS and EMT scores, hypoxia scores, vascular scores and tumor stemness in GSE49710 and GSE62564 datasets, respectively. *** represent p < 0.001.

Figure 7. The correlation analysis between the TLS model and signaling pathways.

(A and B) is the score difference between high and low rating groups in multiple signaling pathways. (C) is enrichment analysis of GO and KEGG. Blue columns correspond to high TLS scores and red columns correspond to low TLS scores. *** represent p < 0.001.

Figure 8. Analysis of clinical characteristics of TLS model.

(A–E) represent subgroup analyses of the TLS scoring model in different INSS stage NB samples, NB samples with or without death, samples with or without MYCN amplification, clinical risk differences, and cytotoxic cytokines. **, *** represent P < 0.01, P < 0.001, respectively.

Figure 9. The construction of Norman map.

(A) is the unifactor COX analysis of clinical features. (B) is the Norman map constructed by integrating clinical factors, and (C) is the K-M curve of the clinical features model.

Figure 10. qRT-PCR and Transwell experiment of NB cell lines.

(A–F) Average expression levels of 6 TLS related genes in two NB cell lines. (G–L) The ratio of the expression level of each gene to the internal reference gene (GAPDH). The red column represents SH-SY5Y and the blue column represents SK-N-(BE) 2. M-R shows the invasion and migration results after WT, NC, and siRNA interference with CD200.