Mutually exclusive teams-like patterns of gene regulation characterize phenotypic heterogeneity along the noradrenergic-mesenchymal axis in neuroblastoma

ABSTRACT

Neuroblastoma is the most frequent extracranial pediatric tumor and leads to 15% of all cancer-related deaths in children. Tumor relapse and therapy resistance in neuroblastoma are driven by phenotypic plasticity and heterogeneity between noradrenergic (NOR) and mesenchymal (MES) cell states. Despite the importance of this phenotypic plasticity, the dynamics and molecular patterns associated with these bidirectional cell-state transitions remain relatively poorly understood. Here, we analyze multiple RNA-seq datasets at both bulk and single-cell resolution, to understand the association between NOR- and MES-specific factors. We observed that NOR-specific and MES-specific expression patterns are largely mutually exclusive, exhibiting a “teams-like” behavior among the genes involved, reminiscent of our earlier observations in lung cancer and melanoma. This antagonism between NOR and MES phenotypes was also associated with metabolic reprogramming and with immunotherapy targets PD-L1 and GD2 as well as with experimental perturbations driving the NOR-MES and/or MES-NOR transition. Further, these “teams-like” patterns were seen only among the NOR- and MES-specific genes, but not in housekeeping genes, possibly highlighting a hallmark of network topology enabling cancer cell plasticity.

KEYWORDS:

• Phenotypic heterogeneity
• noradrenergic-mesenchymal transition
• mesenchymal-noradrenergic transition
• GD2
• neuroblastoma

Introduction

Neuroblastoma (NB) is the most frequent pediatric solid tumor that represents 6–10% of all childhood tumors and accounts for approximately 15% of all cancer deaths in children. It almost exclusively occurs in early childhood; the median age for its diagnosis is about 18 months, and approximately 40% of NB patients are younger than 1 y at diagnosis. It is a neural crest-derived malignancy that manifests along the sympathetic nervous system.^1–3 NB is marked by considerable clinical variability, with the disease spectrum varying from spontaneous regression without needing any treatment to a treatment-resistant tumor exhibiting metastatic dissemination.² Treatment usually includes high-dose chemotherapy, surgical resection, radiation therapy and immunotherapy (anti-GD2 monoclonal antibodies), with most targeted therapy-based inhibitors still in clinical trials.^4,5 Despite intensive therapy, NB patients have high mortality rates, and relapsed patients frequently develop treatment resistance. The median survival for high-risk relapsed neuroblastoma is only 11 months, and there are no curative options for relapsed patients, underscoring the urgent need to improve current treatment regimens.⁶

NB exhibits extensive transcriptional heterogeneity, akin to reports in other cancer types.^7–9 A low mutational burden in the typical NB genome highlights the role of non-genetic heterogeneity in enabling therapy resistance.⁶ Over five decades ago, neuroblastoma cells were shown to exhibit two distinct phenotypes in vitro: the N-type (neuroblast) and S-type (substrate-adherent). The N-type grew as focal aggregates with short neurotic processes and attached poorly to substrate, while S-type were largely flattened, attached strongly to substrate, and resembled non-neuronal precursors. Interestingly, both these cell types could interconvert spontaneously and bidirectionally,¹⁰ reminiscent of phenotypic plasticity reported in other cancers.^11,12 Later, a morphological and biochemical intermediate (I-type cells) were proposed; they had the highest tumor-forming abilities, expressed stem cell markers,^10,13 and were postulated to be a possible precursor to both N-type and S-type cells. Extensive analysis of developmental origins of neuroblastoma have drawn parallels between sympathoblasts – the bipotent cells that can generate both mesenchymal and neuronal phenotypes – and I-type cells. During development, sympathoblasts can give rise to both neuronal (post-ganglionic sympathetic neurons, chromaffin cells) and mesenchymal (mesenchymal stem cells, fibroblasts, Schwann cells, fibroblasts) cell types (Figure 1).¹⁴ Further molecular characterization led to N-type cells being referred to as noradrenergic (NOR) phenotype, while S-type being renamed as mesenchymal (MES) phenotype.¹⁰

Figure 1. Neuroblastoma phenotypes and corresponding developmental states. (top) sympathoblast (I-type) serves as a progenitor of many cell types shown in nor (N-type) and MES (S-type) lineage, at different stages of lineage differentiation. (bottom) nor and MES cells can switch back and forth, indicating reversible cell-state transitions.

Compared to NOR cells, MES cells are resistant to standard chemotherapeutic agents used for neuroblastoma patients – cisplatin, doxorubicin, and etoposide, and are enriched upon treatment of the heterogeneous cell line, SK-N-SH (containing both NOR and MES subpopulations) with these drugs.¹⁰ Consistently, they are also enriched in patients with relapse, as identified by single-cell analysis of primary and relapsed tumors.⁶ This enrichment can be enabled by selection of preexisting MES cells and/or induction of NOR to MES cell-state transition. Single-cell subclones established from NB cell lines could repopulate both cell-states in vitro, indicating reversible dynamic transitions. Such spontaneous stochastic bidirectional interconversion was also seen in vivo.^15–17 Therefore, analogous to the well-studied phenomenon of epithelial-mesenchymal transition (EMT) and its reverse mesenchymal-epithelial transition (MET),^18–20 NB displays phenotypic transitions between the NOR and MES phenotypes – noradrenergic-mesenchymal transition (NMT) and mesenchymal-noradrenergic transition (MNT) (Figure 1).

NMT/MNT can drive the extensive phenotypic plasticity and heterogeneity reported in NB, highlighting the importance of non-genetic mechanisms, such as epigenetic and transcriptomic reprogramming in enabling disease progression and relapse.^16,17,21 While the epithelial-hybrid-mesenchymal spectrum has been thoroughly investigated with respect to its dynamics and impact on hallmarks of cancer,^20,22–24 our understanding of NMT/MNT is relatively limited, and additional research is needed to uncover the underlying patterns of phenotypic plasticity and heterogeneity in NB cell populations. A deeper understanding of the dynamics of cell fate switching in NB could help pinpoint new biomarkers or treatments to target aggressive cell states in NB.

Here, we found a mutually exclusive expression pattern between two sets of transcription factors (TFs) in NB, one driving a NOR phenotype and the other enabling a MES one, thus indicating a “teams-like” behavior between gene expression patterns that determine these phenotypic cell states. This “teams-like” behavior can underlie the phenotypic plasticity and is largely unique to gene lists associated with noradrenergic/mesenchymal axis in NB. We also demonstrated that the NOR- and MES-associated gene lists could reproduce cell-state transitions in NB. Further, our meta-analysis of NB datasets confirmed an antagonistic enrichment of the two phenotypes in bulk transcriptomic datasets along with the enrichment of PD-L1 activity, and glycolytic and fatty acid oxidation (FAO) programs with the mesenchymal state of NB cells. Furthermore, we found that the disialoganglioside GD2, a cell surface marker abundantly expressed in NB tumors, is associated with a NOR phenotype both at bulk and single-cell levels. Our analysis elucidates functional and molecular differences between NOR and MES phenotypes prevalent in NB.

Materials and methods

Software and data and code availability

All computational and statistical analyses have been performed using R (version 4.2.1) and Python (version 3.10.2). Codes used are available at https://github.com/Manas-Sehgal/NB_heterogeneity

Transcriptomic data retrieval and pre-processing

Publicly accessible bulk and single-cell RNA-sequencing datasets were downloaded from NCBI-GEO) repository. Raw gene count matrices were normalized for gene length and transformed to Log2(TPM) (transcripts-per-million) values. Single-cell RNA-sequencing gene-length normalized data (scRNA-seq) was imputed to reduce sparsity using ‘Rmagic’ package.²⁵ Microarray datasets were downloaded from GEO using the GEOquery R Bioconductor package. A total of 75 NB datasets were used for meta-analysis (Table S3). Probe-wise gene expression matrices and their corresponding platform annotation files were used to map the probes to each gene. If two probes matched to the same gene, their mean expression values were used. The resultant gene-wise expression matrices were then log2 transformed.

Principal component analysis and K-means clustering

Principal component analysis (PCA) was used for dimensionality reduction of simulation and transcriptomic data using ‘Scikit-learn’ Python library. K-means clustering was performed on samples in each bulk dataset using the ‘Scikit-learn’ Python library to create two groups (K = 2) based on their segregation along the principal component 1.

For gene swapping analysis, PCA was performed on gene-expression matrices before and after swapping genes from the wild-type 26 NOR/MES gene list was done using ‘prcomp’ function in R. List of housekeeping genes for this analysis was obtained from a previous study (Table S2).²⁶

Gene signatures and gene set enrichment analysis

A total of 14 NOR-specific and 12 MES-specific genes (Table S1) and extended lists of 369 NOR-specific and 485 MES-specific genes from a previous study were used.¹⁷ Gene set enrichment analysis (GSEA) was performed on K-means-generated clusters to examine the enrichment of the noradrenergic and mesenchymal genesets using GSEAPY Python library.²⁷ Gene signatures for Hallmark EMT, fatty acid oxidation (FAO) and glycolysis pathways, Hallmark G2M checkpoint geneset, WP cell cycle, Biocarta cell cycle and KEGG cell cycle were obtained from molecular signatures database (MSigDB).²⁸ Single-sample gene set enrichment analysis (ssGSEA) was performed to obtain sample-wise normalized enrichment scores (NES) for each pathway for bulk transcriptomic data, and AUCell was performed on the imputed scRNA-seq gene matrices.²⁹ Partial EMT and Programmed death-ligand 1 (PD-L1) scores were obtained using previous gene lists.^22,30 All gene lists used in this study are included in Table S2.

Gene correlation matrices and J-metric calculation

Spearman’s correlation coefficients were used to assess the pairwise association between MES-specific and NOR-specific genes. R package ‘ggcorrplot’ was used to plot the correlation matrices. Normalized expression values of genes corresponding to NOR and MES phenotypes (shown in Figure 1) were used to calculate J-metric value for each dataset. To calculate the J-metric value, the sum of Spearman correlation coefficients between genes of different teams was subtracted from the sum of correlation coefficients of genes belonging to the same team. Only significant correlations (p < .05) were considered for this calculation.

Results

NB phenotypic heterogeneity includes the existence of two distinct phenotypes and their corresponding “teams” of genes

To characterize phenotypic heterogeneity in NB, we first collated a set of 14 NOR-specific and 12 MES-specific genes that regulate and/or correspond to noradrenergic and mesenchymal cell-states (Table S1) from a previous study.¹⁷ In our prior studies of epithelial-mesenchymal plasticity in cancer, we have observed that sub-networks of specific genes exist as coordinated “teams” to regulate the dynamics of cellular plasticity.³¹ Players within a “team” activated each other, while those across “teams” inhibited each other directly or indirectly, thus the presence of “teams” can lead to mutually exclusive expression patterns in transcriptomic data. For example, we have pinpointed the epithelial “team” comprised of players such as microRNA-200, OVOL2, GRHL2, and the mesenchymal “team” comprised of players such as ZEB1, SNAIL1 and TWIST.³² Similar “teams” behavior was seen in both the underlying network structure as well as in transcriptomic datasets in additional cancer types, where two mutually antagonistic “teams” promoted invasive/proliferative heterogeneity in melanoma and neuroendocrine/non-neuroendocrine heterogeneity in small cell lung cancer.^7,12 With these features in mind, we sought to understand if the genes involved in regulating the NOR and MES phenotypes in NB were organized in “teams” as well.

To do this, we then computed the pairwise correlations between these genes in multiple bulk transcriptomic datasets (denoted by GSE IDs) comprising NB cell lines and primary tumor samples – GSE9169 (n = 86),³³ GSE17714 (n = 22),³⁴ GSE28019 (n = 24), GSE64000 (n = 8),³⁵ GSE66586 (n = 10),³⁶ GSE78061 (n = 29), GSE19274 (patient samples only; n = 100),³⁷ GSE44537 (n = 6)³⁸ and GSE73292 (n = 6).³⁹ We observed that expression levels of 14 NOR-specific genes were mostly positively correlated with each other but negatively correlated 12 MES-specific genes. Similarly, the 12 MES-specific genes correlated positively with each other, but negatively with 14 NOR-specific ones (Figure 2). This consistent pattern across datasets indicates the existence of two “teams” of players wherein one set of genes corresponds to a NOR phenotype and the other set corresponds to the MES phenotype. Thus, such mutual exclusivity between the two “teams” characterizes NB heterogeneity patterns.

Figure 2. Prevalence of “teams-like” behavior corresponding to the two major phenotypes —noradrenergic and mesenchymal. Correlation matrices illustrating pairwise Spearman’s correlation of NOR-specific and MES-specific genes for GSE9169 (n = 86), GSE17714 (n = 22), GSE28019 (n =24), GSE64000 (n = 8), GSE66586 (n = 10), GSE78061 (n = 29), GSE19274 (n = 100), GSE44537 (n = 6) and GSE73292 (n = 6). All pairwise correlations are significant (p<.05) except the ones indicated by ‘X’. The color bar legend depicts the values of the correlation coefficient. ‘n’ denotes the number of samples in each dataset.

To better understand this “teams-like” behavior, we applied dimensionality reduction to the transcriptomic datasets. In our previous analysis of cell-fate decision networks of EMT and pluripotency, we observed that the underlying “teams-like” structure in a network can reduce the dimensionality of the system and could segregate phenotypes distinctly along principal component 1 (PC1) in a principal component analysis (PCA).⁴⁰ When projected on their first two principal components (PCs), we observed the segregation of samples into two distinct clusters along PC1 across the datasets shown, further highlighting the distinct, “teams-like” behavior of these regulatory pathways (Figure 3a–e). Next, we applied K-means clustering (K = 2) to samples in each dataset to identify two clusters (represented by green and purple color) and calculated the enrichment of the reported 485-gene MES and 359-gene NOR signatures in these two segregated clusters using GSEA.²⁷ We noticed a striking pattern across the datasets, albeit to varying extents, that while one of the two clusters showed enrichment of NOR gene list, the other cluster showed enrichment of the MES gene set (Figure 3a–e, S1). In addition, the extent of enrichment of NOR and MES phenotypes was anti-correlated across datasets (Figure 3f, Table S4). Together, these observations indicate that the two observed clusters in each dataset can be mapped onto the NOR and MES phenotypes.

Figure 3. Two distinct phenotypes—mesenchymal and adrenergic – are observed in NB. a-e) PCA plot showing two distinct classes of samples exist in multiple datasets (left). K-means clustering for K = 2 yields two clusters (shown in green and violet color) along principal component 1 (PC1). GSEA for the mesenchymal (MES) geneset (middle) and adrenergic (NOR) geneset (right) confirms that these two clusters correspond to respective phenotypes for a) GSE9169 (n = 86); b) GSE28019 (n =24); c) GSE64000 (n = 8); d) GSE66586 (n = 10) and e) GSE78061 (n = 29). ‘n’ stands for number of samples in each dataset. Percentage variance explained by each PC is indicated along the respective axes. f). Scatter plot depicting nor (x-axis) and MES (y-axis) ssGSEA scores of samples of GSE66586 (left). Correlation matrix for GSE66586 representing the calculation of J-metric using the correlations between genes shown in Figure 1 (middle) and scatter plot depicting the Spearman correlation coefficient ‘R’ between nor and MES scores across 75 bulk datasets (x-axis) and corresponding J-metric (y-axis) (right). Only 42 datasets in which NOR and MES show significant correlation (p <.05) have been included (Table S6).

Further, we evaluated whether the segregation of samples based on a relatively small gene list (14 NOR-specific and 12 MES-specific) correlated with the clustering obtained through a much larger dimensionality of the system investigated (359-gene NOR signature, 485-gene MES signature). To this end, we used a set of 75 transcriptomic NB datasets for this analysis (Table S3). Previously, we defined a J-metric to quantify the extent of “teams” behavior in each transcriptomic dataset by analyzing the values of a corresponding pairwise correlation matrix.¹² The larger the J-metric, the stronger the “teams” behavior in terms of mutual exclusivity and antagonism between the two phenotypes. Thus, for each transcriptomic dataset, we quantified the J-metric for the pairwise correlation matrix based on 14 NOR-specific and 12 MES-specific genes and calculated the correlation coefficient between the single-sample gene set enrichment analysis (ssGSEA) scores for the 359-gene NOR signature vs. those for 485-gene MES signature. We observed across datasets that the value of J-metric was positive, and the correlation coefficient between NOR and MES signatures was negative. We observed that the J-metric and correlation coefficients between the larger NOR/MES gene lists are negatively correlated with each other (r = −0.33, p < .05, Figure 3f). This trend indicates that both the sets of genes can consistently capture the antagonism between the two phenotypes with a high degree of consistency, thereby suggesting that the “teams-like” behavior also hold true for larger gene-sets belonging to the given biological programs. Put together, these results indicate that two distinct phenotypes can be observed across NB cell lines and primary tumor samples, along with the prevalence of two “teams” of genes driving each phenotype.

“Teams-like” behavior is exclusive to the set of NOR/MES specific genes

Our previous analysis revealed that the “teams” structure was largely unique to the underlying biological network. When compared to an ensemble of random networks of similar size and density and same in-degree and out-degree of each node in the network, the biological network usually had a much higher team strength, thus highlighting that the “teams” topology is a salient feature of networks controlling cancer cell plasticity.^7,12,31 Thus, we next investigated how exclusive the “teams-like” behavior is to the 14 NOR-specific and 12 MES-specific data sets.

As a control, we generated 1000 random ensembles of 26 genes from a set of housekeeping genes (Table S2) and performed PCA for the expression value matrix for those 26 genes in each dataset. We noticed that across datasets, the variance explained by PC1 was much higher for the ‘wild-type’ 26 gene set (14 NOR-specific and 12 MES-specific) than that explained by PC1 corresponding to the 1000 random combinations of 26 housekeeping genes (Figure 4a, i-iii, S4). As another control case, we generated 1000 random combinations of 26 genes with 14 of them randomly chosen from the 369 NOR-specific gene set and 12 of them randomly chosen from the 485 MES-specific gene set and quantified the variance along PC1 for all of them. On average, the ensemble of these 1000 combinations accounted for a higher PC1 variance than the ensemble of 1000 combinations comprising housekeeping genes. This trend suggests that the “teams-like” behavior is exclusive to the specific sets of genes associated with a NOR and MES phenotype.

Figure 4. “Teams-like” behavior is exclusive to the set of NOR/MES team genes. a) i) Histogram (in red) showing the percentage of variance explained by PC1 for 1000 unique combinations of randomly chosen housekeeping genes for GSE28019. The vertical red line depicts the variance explained by the original NOR-MES 26 gene list. Histogram (in blue) showing the same percentage of variance explained by PC1 for 1000 unique combinations based on 14 genes chosen randomly from 369 NOR-specific signature, and 12 genes chosen randomly from 485 MES-specific signature. ii) and iii) are same as i) but for GSE9169 and GSE66586 respectively; b) i) boxplots illustrating distribution of the fraction of variance explained by PC1 (y-axis) vs. number of genes from NOR/MES gene list swapped with housekeeping genes for (x-axis) for GSE28019. The equations for a linear fit of mean values for each swap are also shown. ii) and iii) same as i) but for GSE9169 and GSE66586 respectively; c) heatmap depicting the percentile of variance explained by PC1 by original NOR-MES list in distribution of PC1 variance for random combinations of housekeeping genes; R-squared value (goodness of fit) and Pearson correlation coefficient for linear fit on mean values of variance explained by PC1.

To better understand this trend, we began with replacing one gene at a time in the 26-gene set (14 NOR-specific and 12 MES-specific) with a housekeeping gene. Given many possible combinations depending on which of the 26 genes is (are) being replaced and with which housekeeping gene(s), we tried 100 combinations for each scenario of a fixed number of swaps and plotted the distribution of PC1 variance corresponding to them. We noticed a linear decrease in the average PC1 variance for a given number of swaps, thus, the higher the no. of genes swapped from the 26-gene set, the lower the PC1 variance explained (Figure 4b, S2). This analysis endorses the low dimensionality of the NOR-MES phenotypic landscape and validated the presence of strong “teams” of players exclusive to the genes considered in the analysis. These trends break down once random housekeeping genes that are not connected with NOR-MES axis are introduced.

The trends of 1) a higher PC1 variance of the 26-gene set than a corresponding set of 26 housekeeping genes and 2) a negative correlation between the number of genes swapped and the PC1 variance were both consistently observed across all nine datasets (Figure 4c) investigated (Figure 1).

Experimental perturbations cause phenotypic alterations along the NOR-MES axis

To further substantiate the antagonistic enrichment of the two phenotypes in samples at a bulk level, we obtained ssGSEA scores of each sample across multiple datasets (denoted by GSE IDs) and compared scores of the samples with perturbations to that of the control/wildtype (Figure 5).

Figure 5. Experimental perturbations reveal changes in noradrenergic and mesenchymal programs in bulk RNA-seq data. Barplots depicting experimentally observed significant changes in ssGSEA scores of nor and MES gene signatures in control samples vs. samples with different perturbations in bulk RNA-seq datasets denoted by respective GSE IDs. Error bars represent standard errors in mean values of the replicates and statistically significant differences are indicated by asterisks (*, **, ***, **** for p <.05, <.01, <.001 and <.0001 respectively) for Bonferroni-adjusted p-values.

Topoisomerase-2B (TOP2B) was shown to be essential for SH-SY5Y cells to maintain an adrenergic neural-like (NOR) phenotype, and its deletion upregulated mesenchymal (MES) markers.⁴¹ Consistent with this observation, RNA-seq data for TOP2B-null samples, relative to control, showed a significant reduction in enrichment of a 369 NOR-based gene list and a concomitant enrichment of a 485 MES-based gene list (Figure 5a, i) (GSE142383). Further, ALK inhibitor (TAE684)-resistant SH-SY5Y cells showed AXL activation and induction of EMT.³⁹ Similar behavior was reflected in transcriptomic data where resistant SH-SY5Y cells were enriched in the MES gene set and downregulated in the NOR gene set, when compared to the parental cell-line (Figure 5a, ii) (GSE73292). Next, SK-N-AS cells, upon treatment with the EZH2 inhibitor, tazemetostat, exhibit an increase in NOR scores and reduced MES scores (GSE180512 (Figure 5a, iii). PRRX1A is a MES-specific transcription factor whose overexpression in adrenergic-type SK-N-BE(2C) cells could reprogram them to a mesenchymal state (GSE90804).^17,42 This role of PRRX1A was recapitulated in ssGSEA scores of NOR and MES gene sets (Figure 5a, iv). Along similar lines, knockdown of TCF4, a crucial neurodevelopmental TF, in SH-SY5Y cells led to increased expression of many EMT regulators such as SNAI2, ZEB2 and other TGF-β targets.⁴³ These experimental observations are reinforced in transcriptomic analysis of enrichment of NOR and MES specific gene lists (Figure 5b, i) (GSE48367). Similar results were observed for BE2C cells treated with both JQ1 (a BET bromodomain inhibitor) and THZ1 (a CDK7 inhibitor) that decreased the expression levels of relevant transcription factors, including HAND2, PHOX2B, TBX2, ISL1, MYCN and GATA3.⁴⁴ These six factors are a part of the 14 NOR specific gene set, thus their downregulation via JQ1 and THZ1, as expected, led to a reduced NOR signature and an increased MES signature (Figure 5b, ii) (GSE108914). Further, overexpression of chromatin assembly factor 1 subunit p150 (CHAF1A) – a known inhibitor of neuronal differentiation⁴⁵ – was able to significantly drive up the enrichment of the mesenchymal gene set in SHEP cells after 96 h as compared to control (GSE144311) (Figure 5b, iii). We observed a similar trend with SH-SY5Y cells treated with thapsigargin, an endoplasmic reticulum Ca(2+)-ATPase inhibitor that is known to inhibit cell proliferation.⁴⁶ Given the association of NOR with a higher proliferative status than MES states (Figure S3, Table S5), we noticed a decrease in NOR scores and simultaneous increase in MES scores upon thapsigargin treatment (Figure 5b, iv) (GSE24500). Together, these results highlight how NOR and MES gene lists could capture the cell-state transitions driven by various cellular perturbations.

We next sought to identify additional factors associated with NMT/MNT. First, knockdown of ATRX in NGP (a neuroblastoma cell line) was found to enrich for NOR signature along with a decrease in MES signature (GSE183648) (Figure 5c, i), thus suggesting ATRX as a MES-specific factor. Second, knock out of Huntington’s gene (Htt) in SH-SY5Y NB cells upregulated the NOR gene signature with concurrent reduction in activity of MES signature, as compared to parental cells (Figure 5c, ii) (GSE178467).⁴⁷

Overall, these results indicate that both known and novel experimental perturbations in NB cell lines can drive cell-state transitions along noradrenergic/mesenchymal spectrum in an antagonistic manner.

Meta-analysis of transcriptomic datasets reveals association of NOR/MES phenotypes with metabolic programs and immune markers

Next, we conducted a meta-analysis of 75 bulk transcriptomic datasets (Table S3) using ssGSEA to 1) characterize the modifications in metabolic programs brought on by these cell state transitions and 2) analyze their associations with novel immunotherapeutic targets for NB patients.

First, we noticed that out of 42 datasets that showed significant correlation between NOR and MES signatures, 36 of them (85.7%) were negatively correlated, in support of our previous findings (Figure 6a, left). We next investigated the correlation of NOR and MES genes with a signature associated with PD-L1 activity in carcinomas.²² Consistent with the results observed for carcinomas, out of 37 datasets in which NOR and PD-L1 scores were significantly correlated, 36 (97.29%) showed a negative correlation and 92.3% (36/39) datasets showed a positive correlation with the MES phenotype (Figure 6a, middle, right).

Figure 6. Modes of association between NOR and MES phenotypes with metabolic programs, PD-L1 and GD2 synthase gene (B4GALNT1). a) volcano plots showing Spearman correlation coefficients (x-axis) and -log10(p-value) (y-axis) for NOR vs. MES (left), NOR vs. PD-L1 (middle) and MES vs. PD-L1 gene signature (right). Significant correlations (R > ± 0.3 and p < .05) are shown as red (positive) and blue (negative) datapoints. Same as a) but for b) nor vs. FAO (top), MES vs. FAO (bottom), c) NOR vs. Glycolysis (top), MES vs. Glycolysis (bottom) and d) NOR vs. B4GALNT1 (top), MES vs. B4GALNT1 (bottom). e) scatterplots of scRNA-seq data depicting NOR (x-axis) and MES scores (y-axis) of each cell for the control sample of GSE163429. The vertical and horizontal green lines are positioned at the median of NOR and MES scores, respectively. Maroon datapoints correspond to ‘high’, and blue datapoints represent ‘low’ glycolysis (top), FAO scores (middle), imputed B4GALNT1 gene expression values (bottom). Threshold for ‘high’ and ‘low’ is set at the median value.

Additionally, we evaluated how these two phenotypes associated with two key metabolic processes – glycolysis and fatty acid oxidation (FAO) – that have been shown to be perturbed in NB cells.^48,49 In this context, we observed that in 75.9% datasets, the NOR phenotype is negatively linked with FAO, while the MES signature is strongly positively coupled with FAO (Figure 6b). Similarly, glycolysis scores are also negatively linked with NOR, but show a strong positive link with the MES phenotype. Further, activity of the NOR geneset correlated negatively with the Hallmark EMT signature, while MES scores correlated positively, suggesting similarity between NMT and EMT at a transcriptomic level (Figure S2a). PD-L1 levels were found to be positively associated with both the hallmark EMT and partial EMT signatures, as well as glycolysis and FAO levels across datasets, reminiscent of observations made in carcinomas^22,50 (Figure S2).

Finally, we also analyzed the links between NOR/MES scores with the cell surface marker protein – Disialoganglioside (GD2), which is abundantly expressed on cell surfaces of NB cells and is reported to be a potential target for immunotherapy in NB patients.⁵¹ However, recent observations have indicated that the mesenchymal state of NB cells can confer resistance to anti-GD2 therapy.⁵² To characterize the association of GD2 with NMT, we used the gene expression values of B4GALNT1 (GD2 synthase) – which showed a predominantly positive correlation with the NOR signature and a largely negative correlation with the MES signature, indicating a downregulation of GD2 in cells undergoing NMT. This observation may explain the observation of an enriched mesenchymal tumor population as a hallmark of anti-GD2 therapy resistance in NB patients.

To confirm our observations from bulk transcriptomic data, we also analyzed recently published scRNA-seq data of the NB cell line IMR-575.⁵³ When projected onto the NOR-MES plane, cells showed a significant negative association (R = −0.58, p < 10⁻¹⁵) between the two phenotypes (Figure 6e). These datapoints were then segregated into two groups – high and low – based on their ranks relative to the median AUCell enrichment score for a gene set. We observed that glycolysis^high cells were concentrated at the MES^highNOR^low end of the plane, whereas glycolysis^low cells were predominantly present in the MES^lowNOR^high region (Figure 6e, top). FAO showed a similar trend – FAO^high cells were mostly clustered around the MES^highNOR^low area and FAO^low cells were spread in MES^lowNOR^high portion of the 2D plane (Figure 6e, middle). Furthermore, B4GALNT1 expression was high for cells with MES^lowNOR^high whereas it was low in cells with MES^highNOR^low phenotype (Figure 6e, bottom). These observations are largely concordant with our analysis of datasets at the bulk level and highlight the association of NOR/MES phenotypes with glycolysis, FAO and GD2 expression. These functional synergies might be leveraged in the future to design efficient combinatorial treatment methods that target multiple axes simultaneously to overcome therapy resistance.

Discussion

NB is the most common form of extracranial cancer in infants.¹ Its origin is attributed to improper differentiation of neural crest cells, and the degree of differentiation of these neural crest cells positively correlates with a favorable outcome for patients.^53–56 NB cells are heterogeneous in terms of biochemical and morphological traits, broadly classified into N-type (NOR), S-type (MES) and an intermediate (I-type) cell possessing traits of both N- and S-type cells.⁵⁷ These phenotypes have diverse transcriptomic and epigenetic profiles. While their interconversion has serious clinical implications, the dynamics of NMT/MNT is yet to be thoroughly studied,¹⁰ unlike EMT/MET where systems-level computational and experimental approaches have mapped distinct cell-states, trajectories, and reversibility patterns.²⁴

To unravel the complexity and heterogeneity displayed by neuroblastoma, we collated a cohort of factors that are associated with the two major phenotypes (NOR/MES type).¹⁷ Analyzing their expression levels in multiple bulk transcriptomic datasets pinpointed the presence of mutually exclusive “teams-like” behavior among players that correspond to NOR and MES phenotypes. This behavior is reminiscent of our previous observations in other examples of cancer cell-state plasticity – melanoma, neuroendocrine differentiation and EMT.^7,12,32 PCA led to clear segregation of two clusters, which can be mapped to either a NOR or MES phenotype using GSEA. Importantly, the “teams-like” behavior and consequent PC1 variance are specific to NOR and MES gene lists and were not witnessed for housekeeping genes, indicating functional implications of this mutual antagonism. Such “teams” behavior for epithelial (E) and mesenchymal (M) factors, but not for hybrid E/M ones, was proposed to underlie the difference in plasticity of these phenotypes.³¹ Similar experimental observations for I-type as reported for hybrid E/M phenotypes – higher plasticity, tumor-forming ability and malignancy, and correlation with worse survival as compared to N-type or S-type^13,20,58,59 – and similar “teams-like” behavior seen in transcriptomic data for both EMT/MET and NMT/MNT highlight possible common fundamental design principles in regulatory networks and corresponding phenotypic space. It remains to be seen whether “teams-like” behavior is seen for the I-type as well or if the I-type state emerges because of weakening of “teams-like” behavior among the NOR and MES networks. Future work on identifying I-type specific signatures and comparing I-type, MES-type and NOR-type samples from the same experimental dataset can help us better appreciate the association between “teams-like” behavior and I-type state in NB.

The NOR/MES antagonism was also distinctly observed in our meta-analysis of 75 bulk RNA-seq datasets, and in 10 datasets where we evaluated different experimental perturbations capable of inducing a phenotypic switch along the noradrenergic/mesenchymal axis. Functionally speaking, we observed a positive association of MES phenotype with both glycolysis and PD-L1 signatures, reminiscent of mesenchymal state traits seen in multiple carcinomas.⁵⁰ Expression levels of GD2, a cell surface marker reported to be a therapeutic target for NB patients,⁵¹ correlated negatively with MES, but positively with a NOR phenotype. These trends were also largely consistent at a single-cell level and suggest that combinatorial targeting of GD2 and PD-L1 may be a potent therapeutic strategy to overcome NB heterogeneity. By elucidating the “teams-like” behavior, our analysis indicates how the expression of a set of genes corresponding to a given phenotype is tightly co-regulated, possibly offering additional biomarkers to identify NB phenotypes. Further, we lay the groundwork for future mechanistic mathematical models that can quantitatively explain clinically observed dynamic fluctuations in PDL1 and GD2 expression levels during treatment, potentially due to therapy-driven adaptive plasticity response.

Author contributions

MKJ designed and supervised research. MKJ and JAS obtained funding. MS, SPN and SS performed research. All authors contributed to data analysis and writing of the manuscript.

Acknowledgments

This work was supported by Ramanujan Fellowship (SB/S2/RJN-049/2018) awarded to MKJ by Science and Engineering Research Board, Government of India, and by NCI R01 CA233585 awarded to JAS

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplemental material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15384047.2024.2301802

Additional information

Funding

The work was supported by the Science and Engineering Research Board [SB/S2/RJN-049/2018] and JAS is supported by NCI R01 CA233585.

Notes on contributors

Manas Sehgal

Manas Sehgal graduated from Dr. D.Y. Patil University in 2022 with a Bachelor’s degree in Biotechnology. After a short research stint at Indian Institute of Science, he is now pursuing a master’s degree in Molecular Biosciences at DKFZ in Germany, with a focus on the role of EMT in cancer.

Sonali Priyadarshini Nayak

Sonali Priyadarshini Nayak graduated from the University of Hyderabad with a Bachelors in Systems Biology. She is presently pursuing the Matter to Life Masters program at Max Planck Institute, Germany. She is interested in theoretical biophysics and mathematical modelling of biological systems.

Sarthak Sahoo

Sarthak Sahoo earned his BS & MS in Biology from Indian Institute of Science Bangalore in 2021. Currently, he is a PhD student at the Department of Bioengineering, Indian Institute of Science, Bangalore. His primary research interest is to study the role of phenotypic plasticity and heterogeneity in drug resistance across cancer types using mathematical models.

Jason A Somarelli

Jason A Somarelli is an Assistant Professor at Duke Medical Center and Director of Research for the Duke Comparative Oncology Group. He has over 10 years of experience working in quantifying cellular plasticity in cancers, and focuses on cancer dynamics from a comparative and evolutionary lens.

Mohit Kumar Jolly

Mohit Kumar Jolly is an Associate Professor in Bioengineering at Indian Institute of Science. He has over 10 years of experience in developing mathematical models and conducting transcriptomic data analysis to understand the dynamics of cellular plasticity across multiple cancer types.