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

Data preprocessing of multicenter cohorts

Comprehensive neuroblastoma data were meticulously assembled from two multicenter cohorts, leveraging datasets available in the Gene Expression Omnibus (GEO, accessible at http://www.ncbi.nlm.nih.gov/geo), with specific reference to GSE49710 and GSE62564. Additionally, single-cell sequencing data from neuroblastoma tissues were included from the GSE147766 dataset. These rich datasets provide a robust foundation for our analysis.

Univariate Cox regression

Univariate Cox regression analysis performed using R software (version 4.2.1) with the “survival” package was applied to two separate neuroblastoma datasets (GSE49710 and GSE62564) in our study. Forest plots were generated using the “ggplot2” package (version 3.3.6) for result visualization.

Establishment of a consensus machine learning-driven prognostic signature

To evaluate the consensus machine learning-driven signature (CMLS) and its association with patient prognosis, GSE49710, which contains detailed clinical information, was used as the training set, while GSE62564 served as the validation set. This approach enabled robust assessment of the CMLS’s predictive accuracy across different patient groups. We integrated insights from ten machine learning algorithms, including CoxBoost, stepwise Cox regression, Least Absolute Shrinkage and Selection Operator (LASSO), Ridge regression, Elastic Net, survival support vector machines, generalized boosted regression models, and supervised principal component analysis. This multidisciplinary strategy aimed to enhance the predictive validity and adaptability of the CMLS.

Model optimization for CoxBoost began with the “optimCoxBoostPenalty” function to determine the optimal penalty value, followed by 10-fold cross-validation to identify the ideal number of boosting steps. The “CoxBoost” function was then used for model fitting. Stepwise Cox regression analyses were conducted using the “survival” package, with model complexity assessed via the Akaike Information Criterion (AIC). LASSO, Ridge, and Elastic Net models were implemented using the “glmnet” package, with the regularization parameter lambda determined through 10-fold cross-validation. The alpha parameter, ranging from 0 to 1, adjusted the model between LASSO (alpha = 1) and Ridge (alpha = 0) regularization. The Survival-SVM model was instantiated via the “survivalsvm” function, while the GBM model fitting utilized the “gbm” function, also optimized through 10-fold cross-validation. The SuperPC model was executed using the superpc package, with cross-validation through the “superpc.cv” function. For the plsRcox model, we employed the “cv.plsRcox” function, and the Random Survival Forest model was delineated using the “rfsrc” function, setting “ntree” to 1,000 and “nodesize” to 5.

Prognostic value of CMLS and potential clinical application

Scores were assigned to each sample in both training and validation cohorts based on model outputs, categorizing them into high- and low-CMLS groups. The prognostic relevance of the CMLS was assessed using Kaplan-Meier survival curves, and its prognostic accuracy was evaluated through time-dependent receiver operating characteristic (ROC) curves.

Immune infiltration analysis

In this study, we conducted a comprehensive immune infiltration analysis on 498 neuroblastoma RNA sequencing samples derived from the GSE49710 dataset. To achieve this, we utilized several R packages that are specifically designed for the estimation and visualization of tumor-infiltrating immune cells in gene expression data.

EMT, hypoxia status, tumor cell stemness, and angiogenesis score comprehensive analysis

In this study, we performed a thorough examination of epithelial-mesenchymal transition (EMT), hypoxia, tumor stemness, and angiogenesis in tumor samples. We assessed EMT using 200 genes from the HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION pathway and calculated EMT scores for each sample via single-sample Gene Set Enrichment Analysis (ssGSEA). Hypoxia was evaluated with ssGSEA on genes from the HALLMARK_HYPOXIA pathway to assign hypoxia scores. We broadened our analysis to include KEGG-listed pathways to understand the activity of various signaling pathways within tumors, calculating their activity scores with ssGSEA. For tumor cell stemness quantification, we utilized human stem cell datasets from the Progenitor Cell Biology Consortium (PCBC) on Synapse.org and applied the One Class Linear Regression (OCLR) method to measure stemness levels.

Real-time quantitative PCR and transwell experiment

NB cell line was purchased from Meisen Cell Biotechnology Co., LTD., Zhejiang, China. For the primer sequences corresponding to each gene in qPCR experiments, detailed information can be found in Table S1 (submitted). Both the SK-N-BE (2) and SH-SY5Y neuroblastoma cell lines were amplified using standard cell culture techniques. In order to compare the relative expression levels of the two cell lines, GAPDH was used as the internal reference gene, and the other six genes were used as the experimental group. First, the culture solution of the petri dish was sucked out, and then the cells were cleaned once with 1 × PBS and the adherent NB cells were removed with a cell scraper. After adding 1 mL of RNAex and shaking well, all liquids were transferred to an enzyme-free centrifuge tube for centrifugation, and placed in an ice box for 5 min of precipitation to prepare for subsequent purification. Before extraction and purification, 1 ml cold ethanol of equal volume (70% concentration) was added, mixed, transferred to RNAmini column and centrifuged at 12,000 rpm for 1 min. Discard the bottom liquid that has been centrifuged, add 600 ul Buff RW1 and mix well. Then centrifuge under the same conditions for 1 min and remove the bottom liquid. Next add 650 ul Buff RW2 to clean twice, then transfer to the collection tube with column core and centrifuge. After removing the lower liquid, the column core was transferred to a new collection tube, and 40 ul of enzyme-free water was added to the middle and placed on ice for 30 min. Centrifuge at 12,000 rpm for 3 min, and then the RNA concentration was measured. According to the indication of the reverse transcription kit we used, the RNA concentration should not be lower than 70 ng/mL. mRNA was extracted using a commercial kit (Total RNA Purification Kit) and quantified with a spectrophotometer. This kit was purchased from Accurate Biotechnology (Hunan, China) Co., Ltd. Using a NanoDrop spectrophotometer, a blank measure is first used and then the concentration of the sample is tested. The ratio of A260/A280 was calculated, and the ratio was ≥1.8, which met the experimental requirements. According to the instructions, 30 ul system is used for reverse transcription. According to the instructions of the RNA-to-cDNA™ Kit(Applied Biosystems™), firstly, add 6 ul of 5× primescript buffer and 1.5 ul of ENZYME MIX, Olido dT primer and Random 6 mers each. Then add no more than 19.5 ul of RNA samples. Finally, RNA-free H2O was supplemented to make the system reach 30 ul. Reaction conditions: Reverse transcription temperature 42 degrees Celsius for 15 min, inactivation temperature 85 degrees Celsius for 5 min total of 20 min. The 20 ul qPCR system was respectively added with: 10 ul of 2X SYBR Green Pro Taq HS Premix, 2 ul of CDNA, 1 ul each of the forward primer and reverse primer, and 6 ul of RNAse free H2O, according to the reagent kit (Accurate Biology, Guangzhou, China). RT-PCR experiments were conducted on the Bio-Rad CFX Manager platform. The GAPDH gene served as an internal reference gene, with each target gene having three separate reaction wells to ensure the accuracy of the results.

si-CD200 was designed and customized according to the siRNA sequence. siRNA was diluted to an appropriate concentration (20 nM) and mixed with Lipofectamine RNAiMAX transfection reagent, and the transfection complex was added to the culture dish of NB cells and incubated for 48 h. The concentration of si-CD200 is based on the possible optimal concentration range obtained from the previous concentration gradient pre-experiment, which is between 10 and 30 nM. The culture medium was DMEM medium with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin (P/S) mixture. Incubation was carried out in an incubator at 37 °C and 5% CO2. An appropriate amount of basement membrane (Matrigel) is applied to the porous membrane of the Transwell chamber and hydrated to simulate the extracellular matrix. Place the Transwell chamber into a 24-well plate. The cell suspension including 1 × 10⁵ was added to the upper chamber of Transwell chamber and the medium was added to the lower chamber. The Transwell chamber is cultured in a cell incubator to give the cells time to migrate or invade through the porous membrane. After culture, migrating or invading cells were immobilized with paraformaldehyde and then stained with crystal violet stain. By looking under a microscope and counting the cells that migrate to the underside of the membrane.

Single cell sequencing analysis and correlation analysis of drug sensitivity

To begin with, we implemented stringent quality assurance measures on the acquired GSE147766 dataset. This was essential to guarantee the integrity and dependability of the data. In this process, we evaluated the sequencing depth to ensure it met coverage requirements, discarded cells that did not comply with quality criteria, and eliminated cells containing a disproportionately high percentage of mitochondrial genes (Fig. S1). To perform principal component analysis (PCA) on our preprocessed dataset, we utilized the RunPCA function from the Seurat package. We then reduced the dimensionality by selecting high variable genes (HVGs) and focusing on the eighth principal component (PC8). Following this, we employed the Rtsne function to conduct t-distributed stochastic neighbor embedding (t-SNE) clustering analysis. Finally, we annotated the clustered cell populations using CellMarker 2.0 to identify the distinct cellular components. After annotation, the expression of target genes in UMAP is visualized and grouped using violin maps. Finally, we divided 17 NB samples into newly diagnosed group (14 NBs) and refractory relapse group (3 NBs) according to clinical data for single cell level visual analysis.

Furthermore, the relationship between the mRNA levels of the six candidate genes and the sensitivity to common anticancer drugs was analyzed using pan-cancer data from the Genomics of Drug Sensitivity in Cancer (GDSC) and Cancer Therapeutics Response Portal (CTRP). These analyses were performed on the Gene Set Cancer Analysis (GSCA) platform at wchscu.cn.

Statistical methods

We conducted statistical analyses and generated graphs utilizing R software (version 4.3.3). To evaluate the differences between the two groups, we employed both two-tailed unpaired Student’s t-tests and Wilcoxon rank-sum tests. For the comparison of categorical variables, Fisher’s exact test was employed. The threshold for statistical significance was established at P < 0.05. The above statistics were analyzed using SPSS 22.0 software (IBM Corp., Armonk, NY, USA).

Construction of an efficient prognostic model based on tertiary lymphoid structure-related genes in neuroblastoma

We conducted a thorough literature review to identify genes associated with TLS in cancer, yielding 37 TLS-related genes (Sautes-Fridman et al., 2019). This gene set encompasses 12 chemokine signature genes, eight T follicular helper (TFH) cell signature genes, 15 T helper 1 (TH1) cell and B cell signature genes, and two plasma cell and CXCL13 signature genes (Sautes-Fridman et al., 2019). As depicted in Figs. 1A, 1B, we performed univariate COX regression analysis on the expression levels and prognostic data of these 37 TLS-related genes in two neuroblastoma transcriptome datasets, GSE49710 and GSE62564. Genes to the left of the forest plot line with a p-value less than 0.05 were considered to be NB prognostic TLS-related genes.

Establishment of a consensus machine learning-driven prognostic signature

We utilized the GSE49710 dataset as the training set and the GSE62564 dataset as the validation set to incorporate the prognostic TLS-related genes identified through univariate COX regression analysis into a consensus machine learning-driven algorithm comprising ten distinct machine learning methods. The goal was to select the machine learning method with the highest concordance index (C-index) to construct the prognostic model. As shown in Fig. 2A, the random forest model ranked fifth with the highest C-index of 0.955. Analysis of the optimal number of features to include, as recommended by different machine learning combinations, revealed that the random forest model achieved the best accuracy when six feature genes were included (Fig. 2B). The area under the curve (AUC) also indicated that the performance plateaued when six features were incorporated (Fig. 2C). Consequently, we developed a prognostic scoring model based on six TLS-related genes: Risk score = (−0.274148539440277 * expression of CCL4) + (−0.503939241290335 * expression of CD200) + (−0.356039730964334 * expression of CXCR3) + (−0.109685444318236 * expression of CCL2) + (0.195937517824706 * expression of CCL21) + (0.0498379732890246 * expression of IGSF6).

Gene Ontology (GO) analysis revealed that these differentially expressed genes are enriched in multiple pathways involving cytokine and chemokine receptor activation, and they participate in the composition of membrane signal receptor complexes, mediating leukocyte migration and various chemokine responses (Fig. 2D).

Association of TLS-related genes and the constructed model with neuroblastoma patient survival

To validate the efficacy of the feature genes selected by machine learning, we employed the GES49710 dataset as a training set for single-gene efficacy validation. As shown in Figs. 3A–3F, high expression levels of CCL2, CCL4, CCL21, CD200, CXCR3, and IGSF6 were significantly associated with improved survival in neuroblastoma patients (all p < 0.001). These trends were also observed in the independent external validation set GSE62564 (p < 0.001, p < 0.001, p = 0.007, p < 0.001, p = 0.005, and p = 0.002, respectively) (Figs. 3G–3L). Subsequently, we further validated the prognostic power of the TLS gene score model in predicting neuroblastoma event-free survival (EFS) and overall survival (OS). Results in Figs. 4A, 4B, 4G and 4H demonstrated that the low-risk score neuroblastoma group exhibited significantly better EFS and OS in both the training and validation sets (all p < 0.001). Diagnostic ROC curves confirmed the predictive efficacy of the score model for OS and EFS (AUCs of 0.943 and 0.747; 0.795 and 0.681, respectively) (Figs. 4C, 4D, 4I, 4J). Time-dependent ROC curves also showed that the score model had high predictive efficacy for OS and EFS at 0.5 and 1 year (AUCs ranging from 0.611 to 0.979) (Figs. 4E, 4F, 4K, 4L).

Immune infiltration scores

We employed multiple assessment methods, including TIMER, ESTIMATE, and CIBERSORT, to analyze differences in immune cell infiltration between groups with high and low model scores. Figure 5A presents the relative abundance of six immune cell types: B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils, and dendritic cells, indicating that the low model score group has a higher abundance of various immune cells. ESTIMATE results revealed that the low-score neuroblastoma group had higher stromal, immune, and ESTIMATE scores, implying a lower tumor cell purity (Figs. 5B–5D). In contrast, the CIBERSORT algorithm score showed that the low-score group had a higher proportion of Tregs and lower proportions of CD8+ T cells and NK cells. Consistently, the distribution of macrophages and B cells was similar to the other two immune scores, with the low-score group having a higher proportion of high-function immune cells such as M1 macrophages (Fig. 5E). The overall correlation analysis is displayed in Fig. 5F.

Comprehensive scores for EMT, hypoxia, angiogenesis, tumor stemness, and signal transduction pathways

EMT, hypoxia, angiogenesis, and tumor stemness are key indicators for evaluating tumor recurrence, metastasis, and patient prognosis. We calculated these scores using the ssGSEA algorithm, and the results, as shown in Fig. 6, indicate that the neuroblastoma group with lower model scores has higher EMT, hypoxia, and angiogenesis scores (Figs. 6A–6C, 6E–6G) and lower tumor stemness scores (Figs. 6D, 6H). Tumor-related signal transduction pathway scores revealed that the group with lower scores has higher scores for various signal transduction pathways (Figs. 7A, 7B), which is consistent with the results of KEGG and GO pathway enrichment (Fig. 7C).

Clinical feature analysis and nomogram construction

Subsequently, we conducted a detailed analysis of the relationship between the scoring model and clinical characteristics. As shown in Fig. 8A, patients with International Neuroblastoma Staging System (INSS) stages 3 and 4 had higher scores, which is consistent with the higher scores observed in clinically high-risk and deceased patient groups (Figs. 8B, 8D). MYCN amplification is known to predict poor prognosis in neuroblastoma (NB) patients, and our subgroup analysis revealed that the NB group with MYCN amplification had higher scores (Fig. 8C). Given the potential association between tertiary lymphoid structure formation and immune cell function, we analyzed the relationship between model scores and cytotoxic factors. The results, as depicted in Fig. 8E, showed that the low-score NB group had higher levels of GZMB, PRF1, IFNG, IL-2, IL-15, and IL-18. Based on these findings, we included multiple clinical characteristics such as gender, age, staging, and MYCN amplification in univariate COX regression and constructed nomograms suitable for clinical practice that incorporated factors significantly related to NB prognosis (Figs. 9A, 9B). Subgroup analysis of the newly constructed clinical characteristic scoring model revealed that the low-score group had significantly longer survival times compared to the high-score group (Fig. 9C).

qRT-PCR and transwell experiment of NB cell lines

In order to detect the difference in expression of six TLS signature genes in non-MYCN amplified and MYCN amplified cell lines, we selected SH-SY5Y and SK-N-BE (2) for real-time quantitative PCR detection after conventional cell culture amplification. Results as shown in Figs. 10A–10F, after three multiple hole repeats, except for the relative expression levels of CCL21 and IGSF6 increased in SH-SY5Y cell line, the other four characteristic genes were significantly lower than SK-N-BE (2) cell line. Figures 10G–10I shows the difference in the expression of each characteristic gene relative to the mean ct value of GAPDH gene in two different cell lines. Except for IGSF6, the expression level of other characteristic genes in SK-N-BE (2) was higher than that of SH-SY5Y. Considering that CD200 was assigned the largest coefficient in our model and that we found the phenomenon of CD200 aggregation and expression in tumor cells through single-cell cluster analysis, we decided to study the effect of CD200 on NB cells. Next, we down-regulated CD200 expression in NB cells using small interfering RNA, and then verified the effect of CD200 on invasion and migration of NB cells. As shown in Figs. 10M–10R, the invasion and migration ability of both NB cell lines decreased significantly than WT group and NC group after interference with CD200.

Single-cell sequencing analysis of NB samples

To gain a deeper understanding of the expression patterns of TLS model-related genes across different immune cells, we performed tSNE and UMAP clustering analysis and annotation on single-cell sequencing samples from 17 NB tissues (Figs. S3A–S3D). Figure S3E presents a bar chart showing the proportion distribution of various immune cells in the 17 NB samples. Further single-gene UMAP analysis revealed that CCL2 and CCL21 are primarily distributed in smooth muscle cells, tissue stem cells, endothelial cells, and a subset of macrophage subpopulations. CCL4 and CXCR3 are mainly found in NK cells, T cells, and monocyte-macrophage subpopulations. CD200 was expressed in various immune cells to varying degrees. In our single-cell cluster analysis, we found a very interesting phenomenon: CD200 was expressed at a higher level in neuron-like tumor cell subpopulations. IGSF6 is mainly concentrated in monocyte-macrophage subpopulations (Figs. S3F–S3Q). Subsequently, we conducted subgroup comparative analysis on single-cell samples from three refractory relapsed NB cases and 14 newly diagnosed NB cases (Figs. S4, S5). Figures S4, 5A–5D show the cellular distribution maps after tSNE and UMAP clustering analysis and annotation of RR-NB and ND-NB single-cell samples. Figures S4, S5E presents the statistical results of cellular distribution differences in different NB subgroup samples. It is evident that, compared to most ND-NB samples, RR-NB samples have a higher proportion of neuron-like NB cells and a lower proportion of T cells and B cells (Figs. S4, S5E). Single-gene UMAP distribution and statistical analysis found that the enrichment density of the six characteristic genes in the 14 ND-NB sample cell subpopulations was significantly higher than in RR-NB samples (Figs. S4, S5F–S5Q).

Correlation analysis of anti-tumor drug sensitivity

Finally, we utilized the Genomics of Drug Sensitivity in Cancer (GDSC) and Cancer Therapeutics Response Portal (CTRP) databases to perform correlation analysis between single-gene expression and drug sensitivity. The results showed that the transcriptomic expression level of CCL2 is significantly positively correlated with the sensitivity to the TOP30 anti-tumor drugs, while the other five genes exhibit varying degrees of negative correlation with most drugs. Notably, the expression levels of IGSF6, CXCR3, and CCL4 genes showed the most pronounced negative correlation with drug sensitivity (Figs. S2A, S2B).